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THE SIZE EFFECT OF FERRITE SPHERESt+ 
By WALTER HAUSER and LOWELL BROWN 
(Lincoln Laboratory, M.I.T., Lexington, Massachusetts) 
[Received 14 May 1959] 


SUMMARY 

An integral equation technique is used to obtain an approximate formula for the 
apparent susceptibility of a small ferrite sphere of radius a within a microwave 
cavity. The first-order effect of a plane wall is also considered using the quasi- 
stationary approximation that the fields within the sphere are uniform. It is shown 
that for spheres at a distance of the order of a wavelength A away from the cavity 
wall, the wall contributes to the apparent susceptibility of the sphere terms of the 
order of (ka)® and higher, where k = 27/A, justifying to order (ka)® the use of free 
space calculations in cavity problems to estimate the fields within the sphere. 


1. Introduction 

Tue problem of the effect of a small ferrite sphere on the resonant fre- 
quency of a microwave cavity has attracted the attention of several 
investigators (1). The computation of the shift in resonant frequency or 
the apparent susceptibility of the sphere involves a knowledge of the 
electromagnetic field within the perturbing sphere. An exact solution of 
this problem not being possible, we are obliged to seek approximate 
solutions. The most accurate approximate solution has been obtained by 
Hurd (1), who utilized a power series solution in A~! for the scattering of 
a plane electromagnetic wave from a small anisotropic ellipsoid in infinite 
space to approximate the field within a sphere of radius a to order (ka)?, 
where the wave number k = 27/A. The degree of accuracy of the apparent 
susceptibility of a ferrite sphere within a cavity computed on the basis 
of fields within the sphere in infinite space remains however an open 
question. Dnestrovskii (2) in a very interesting way has shown that the 
quasi-stationary approximation leads to results which are good to order 


(ka)?. 


In this paper we utilize an integral equation technique in order to 
obtain an approximate solution to this problem. With the introduction 
of the free space dyadic Green’s function, we are able to obtain a formal 
solution to the problem in terms of integrals involving the field vectors 
in the perturbing sphere plus integrals over the surface of the cavity. 


+ The work reported in this paper was performed by Lincoln Laboratory, a centre for 
research operated by Massachusetts Institute of Technology with the joint support of the 
U.S. Army, Navy, and Air Force. 


(Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 3, 1960) 
5092.51 Ss 
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Addition of the integral equation satisfied by the unperturbed cavity mode 
permits us to show that some of the remaining integrals contribute terms 
of order (ka)® and higher when the sphere is a distance the order of a 
wavelength from the cavity walls. Neglect of these integrals yields the 
integral equation for a sphere in infinite space having the unperturbed 
cavity mode incident upon it. We are thus able to justify to order (ka)’ 
the use of free space calculations such as Hurd’s (1) in cavity problems. 

tecasting the free space integral equation in the form of a variational 
principle, we proceed to obtain a first-order estimate of the effect of a 
plane wall and the effect of the size of the sphere on the apparent suscepti- 
bility of the sphere. 


2. Formulation of the problem 

We consider the problem of an anisotropic obstacle within an ideal 
cavity free of charges and currents. The fields within the cavity satisfy 
the differential equations 

VxE = —jwp,H—jop’.H, 

VxH = jwe, E+jwe’.E, (2.1) 
where ¢, and yu, are respectively the unperturbed scalar electric permit- 
tivity and magnetic permeability. The dyadic electric permittivity € and 
the dyadic magnetic permeability p of the obstacle are given by 

e=e,I+e’, 


p= pol+p’, (2. 
where I is the unit dyadic. e€’ and p’ are to be considered as functions of 
position assuming the value zero outside the obstacle. 

Equations (2.1) together with the boundary conditions that the normal 
component of B and the tangential component of E vanish on the cavity 


wall, n.B=nxE=0 onS, (2.3) 


pose the eigenvalue problem for the resonant frequency, w = kc, of the 
cavity. 


The two equations (2.1) can conveniently be expressed in the single 
matrix vector equation 


[o,V x —k@.JF = 


where rs 
. | h)’ 


Og = (" a and UM = & ‘) 
j 9 0 B/Ho 


In the absence of the perturbing obstacle the solutions to equation (2.4) 
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are the normal modes of the cavity, 


which satisfy the differential equation 
[o,.V x —k, F .|F, = 9, 


where JF = ~ @ 


and k,, are the resonant wave numbers of the cavity. The complex con- 
jugate transposed normal modes, 


F}, = (ef, hy), (2.8) 
satisfy the differential equation 
V x (Fio.)—k, Fi = 0. (2.9) 
We normalize the modes so that 


| F},.F, dr = { (eke, +h¥.h,) dr = 2, (2.10) 


where the integration extends over the volume V of the cavity. 

Insertion of a small obstacle into the cavity produces a change in its 
resonant frequencies as expressed by the Bethe—Schwinger perturbation 
formula (3) which we now proceed to derive. 

Multiplying equation (2.4) from the left by #}, equation (2.9) from the 
right by F, subtracting the results and integrating over the volume of 
the cavity yields 

([(W x Fh)oy. F—F},.0,(V x F)] dr 


== (k,—k) | F}.F dr—k | F).U'.F dr, (2.11) 
where @’ = U—SF. 
Now (V x F})o,. F —F},.0,(V x F) = V.[(F} 0.) x F], 


and utilizing the divergence theorem and the boundary conditions satisfied 
by F and F},, we find that 


[ V.(Fh0,xF) dr = J (nx F!).0,F dS 


= jf [—(axe%).h+(mxht).e]dS = 0. (2.12) 
It follows therefore from equation (2.11) that 
k,—k | F}.'.F dr 
kf FLLF dr 
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To solve equation (2.4) we introduce the matrix dyadic Green's func- 
tion, G(r r’), which satisfies the differential equation 
(o,.V x —kF .)G(rir’) = d(r—r)F, (2.14) 
where 5(r—r’) is the three-dimensional Dirac delta function. Writing 
G(rir’) = (0.V x +kF.)T(rr’), (2.15) 
it follows from equation (2.14) that T(rir’) satisfies the differential 
equation (Vx Vx —k)F(rir’) = &(r—r')F. (2.16) 


Thus we find that T(r r’) is a diagonal matrix dyadic Green’s function 
satisfying the inhomogeneous vector wave equation. If the electric and 
magnetic parts of the matrix dyadic Green’s function satisfy respectively 
the homogeneous boundary conditions imposed upon the electric and 
magnetic field, we find that 


T(r r’) M(rir) és ); 

0 M(r r’) 
where N(rr’) and M(rr’) are respectively the electric and magnetic 
dyadic Green’s functions (4). 

In the work to follow we shall make use of the free space matrix dyadic 
Green’s function, G(r r’), for which (5) 


I(r r’) g.|I EVd lar r’), (2.18) 


where g(r r’) is the scalar Green’s function satisfying the differential 


equation V°g(r ir’) +k¢(r vr’) = —d(r—r’), (2.19) 


as can readily be veritied by substituting equation (2.18) into (2.16). For 
our problem of a ferrite-sphere within a cavity, we choose for g(r r’) the 
standing wave solution 
cosk r—r’ 
g(rr) . 
4nir—r 


whose expansion in spherical coordinates is 


= » l ! , 
g(r ir’) SY aeypys Em. m): cos m($—¢’) {'(cos @) X 
do ~- — (l+-m)! 


| jikrynfkr’) ip = 


< P(cos 6 Lee 
l (cc ) ni kr) j (kr ) (r 


where = 


jl m=0 

(2 m9, 

P*"(cos @) is the associated Legendre polynomial of the first kind, and j,(z 
i a) } . Ji 

and n,(z) are the spherical Bessel and Neumann functions. 
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In addition, for the wall effect we shall require the matrix dyadic Green’s 
function G,,(r r’) satisfying the above-mentioned homogeneous boundary 
conditions on an infinite perfectly conducting plane. The electric and 
magnetic dyadic Green’s functions for this problem have been given by 
Levine and Schwinger (5). We find that 
P(eir’) = Prer’)+Pir'r’), 


where rrr’) = a, avy" 1. -2nn|g(r ir’). 


a=(, 4 


n is the unit normal to the plane, and 
r”°=r’—2(n.r’)n 
if the origin of the coordinate system lies on the plane. 

In dealing with matrices whose elements are dyadics we have to distin- 
guish between the operation which transposes both the matrix and its 
dyadic elements, and the operation which transposes only the dyadic 
elements. The first we shall designate by the superscript 7’, the latter by 
the superscript ¢. For example, if 

ab cd 
FL , 
ef gh 
ba _ fe 
dc hg)’ 
# - dc 
fe hg 
and "= g. 
We find that #7” satisfies the differential equation 
V x Flop —kFT—kF™T.MU'T = 0. 


then GFT — 


We now introduce the adjoint matrix dyadic Green’s function, G?(r r’), 
which satisfies the differential equation 


Vk Grr’ )—k@iir rv’) = d(r—r')F. 2.29) 
For real k Grr’) = @*(r'r’). 
In general, if the proper boundary conditions are imposed upon the matrix 
dyadie Green's function, 
Grr’) = G7 (rr). (2.30) 
To show this, multiply equation (2.29) by the arbitrary matrix vector a 
and the equation satisfied by @(rir”) by the arbitrary vector b. We thus 
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obtain the equations 
o3V x Gi(rir’)—kG@i(rir’) = d(r—r’)a 
and V «Gf(rir’)of—kG@i(r ir’) = b73(r’—r’), 
where Gi(rir’) = G@irir’).a 
and Gi (rjr”) = [G(r\r’).b]?. 
Multiplying equation (2.31) from the left by @7(r\r’), (2.32) from the 
right by @2(r r’), subtracting the results and integrating over the volume 
of interest vields 
b?.[G7(r'\r’)— G(r’ r’)].a 

-| n.GF(r r’)o} x Gi(rir”) dS. (2.33) 
The right-hand side of this equation vanishes if the elements of the matrix 
dyadic Green's function satisfy the proper homogeneous boundary condi- 
tions, thus yielding the reciprocity relation (2.30) satisfied by the matrix 
dyadie Green’s function. 

Now, multiplying equation (2.28) from the right by the adjoint matrix 
dyadic Green's function, (2.29) from the left by #7, subtracting the results 
and integrating over the volume of the cavity yields 

Fir’) =k | Ar). WT. Girir’) dr— | nx FF? (r).ofG(r\r’) dS. 

‘ (2.34) 
Similarly for the cavity mode .¥,, which satisfies equation (2.4) with 
,  k—k 
U = - we 


we find 


(ky k) [ FJ (r) Gr r’) dr— 
( (nx Fe (r)).ofGi(rir’) ds. 2.35) 


Subtracting the transpose of equation (2.34) from the transpose of (2.35) 
and using the reciprocity relation of the Green’s function yields 


Fe’) = Alr')+k | Gr’ rv). WF (vr) dr+ 


+ (Ky—k) ( G(r’ r).A,(r) dr— ( G(r’ r)o,.(nf(r))dS, (2.36) 


where f(r) = ¥(r)— A,r). 

We shall now show that if @(r’'r) is the free space matrix dyadic 
Green's function, the last two integrals in equation (2.36) contribute only 
terms of the order of (ka)? and higher for obstacles placed at a point ro 
near the centre of the cavity. 
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From the Bethe-Schwinger perturbation formula (2.13) it follows that 
as 
(ky—k) = ky x Oo *), 
Furthermore, we find that 
[ Gj(r'\r). A(x) dr 


i 2.0" x +; vv'+ e| If gleo\e) Flr) dr 


approaches k | g(to\r) Ar) dr = Of“) Fama 


for small k. R is a distance of the order of a cavity dimension, which is 
of the order of 1/& for the lowest cavity modes. Thus we see that 


(ky—k) [ Gylroir).A(r) dr = O(ka)*. (2.37) 


As for the last integral in equation (2.36), to a first approximation f(r) 
is the field of the electric and magnetic dipoles induced within the sphere 
by the cavity field A, plus the field of the image dipoles. From the 
results of the next section we shall see that for a sphere a distance d ~ 1/k 
away from a plane wall of the cavity, the contribution of the last integral 
is of the order (ka)3. This is because the lowest order of the induced 
moments, which are of the order AV = O(a’), produce a field 


1 ik 


f(r) = | — |e = O(ka)s 
rij} J 


\ 


ro—r iro—re iro 
t_ _ oh). 


0 
Thus for obstacles at a distance the order of a wavelength away from 
the walls of the cavity 


( Gri \r)oy.(n x f(r)) dS = O(ka)s| [ Aro r) dr| = O(ka)3. 


Hence, for an obstacle located near the centre of the cavity, 
Fr’) = Ale’) +k | G(r'\r).W' .F(r) dr+-O(ka)s. (2.38) 
Neglecting terms of order (ka)® yields the integral equation for an 
obstacle in free space with an incident field A,(r). This justifies to order 
(ka)* the use of the free space calculations such as Hurd’s (1) in cavity 
problems. A similar analysis shows that for the centre of the sphere near 
a plane wall of the cavity but far from other walls 


Fr’) = Fir')+k | G(r’ |r). @' .F (vr) dr+O(ka)3. (2.39) 
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We shall now proceed to obtain solutions to equations (2.38) and (2.39) 
within the sphere in the form of a power series expansion 


F(r)= A+B.r-+..., (2.40) 


where A is a constant matrix vector and B a one-column matrix dyadic, ... . 
The coefficients A, B,... are readily determined by the use of a variational 
principle. Utilizing the normalization condition (2.10), we find that the 
variational expression 


S) ky 


k , ‘ee 
, + | FF... Fi dr 


| FU F dr+k || FM'S e'\v).U'.F drdr’ (2.41) 


vields integral equations (2.38) and (2.39) upon first-order variation with 
respect to the fields F*. 


3. Wall and size effect for ferrite spheres 

A. Wall effect 

As a first application of equation (2.41) we consider the case of a very 
small sphere (ka << 1) for which we can approximate ¥(r) within the 
sphere by the first term of equation (2.40), the constant vector A, Substi- 
tuting the approximate fields 


¥(r) =A, Fir) = At 


into variational expression (2.41) and varying with respect to A’ yields 
the integral equation 


A = F,(ry)4 “7 | | G(r’ r).W’. Adrdr’ + O(ka)?, (3.1) 


where the integration extends over the volume of the sphere AV, and ry 
is the position vector of the centre of the sphere. 

If the sphere is far from the cavity walls we use the free space Green’s 
function for G(r’ r). From equations (2.15) and (2.18) we find that 


Grr) — |o,' xs vv Eg | ae r), (3.2) 


and using the expansion of g(r’ r) in spherical coordinates (equation 
(2.21)) we obtain 


— | | G(r vr) drdr’ [1 +-2(ak)n,(ka)j,(ka)|F ~ —)F + O(ka)?. 
7 (3.3) 
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Hence, to order (ka)?, 
A = [F+4%')".F,(r,). (3.4) 


Substituting this approximation for the fields within the sphere into the 
Bethe-Schwinger formula (2.13) yields 


i Ft). MU! [IFAW Fly). (3.5) 
This approximation is referred to as the quasi-stationary approximation 
and is seen from our analysis to be valid to the order (ka). This estimate 
of the order of accuracy of the quasi-stationary approximation has also 
been obtained by Dnestrovskii (2). 

If the sphere is located at a small distance d from a plane wall of the 
cavity, but far from other walls, we must use the Green’s function G,,(r’ \r) 
defined by equations (2.15) and (2.22). Assuming once again uniform 
fields within the sphere we must now evaluate 


2 | G(r’ ir) drdr’ . | [ (60 r)+@G,(r'\r)| drdr’. 
AV J. AV J, 
The first integral on the right is given by equation (3.3) and the second is 
quite simple to evaluate. Within the region of integration @,(r’ |r) has 
no singularities and thus satisfies the scalar Helmholtz equation both as 
a function of r and r’. Now any function V(r) which satisfies the scalar 
Helmholtz equation has an expansion in spherical coordinates given by 
the equation 


V(r) . D2 Ay MAT)Y (8, ¢). 


" 


Thus 


[ V(r) dr = drrdqq [ jglkr)r? dr = 4rdq9j,(ka)a® k 
1 0 


3 
=< ool 1 —j,(ka)? + O(ka)'}. 


But Ao = ¥ (9), 


therefore, to order (ka)', 
3 
{¥) dr = = [1— (ka)? }¥'(0). (3.6) 
Applying this mean value theorem to G(r’ |r) yields, to order (ka)?, 


< | | G(r’\r) drdr’ = (AV)KG,(r4, |r) 
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which from equations (2.15) through (2.25), to order (ka)?, becomes 


Pe 1 /a\3 é — 
n>» 1+5(5) o,I+nn], (3.7) 


wire oe Es (ka){a 
ar | | ner eer eas (; 


(Oo WV aa 
JP = 0293 = i(, a) (3.8) 


Inserting (3.3) and (3.7) into equation (3.1) we obtain 


; : \3)-1 
A = [94.39 42% (kay(2) n <M —* 01mm). (4 F(t), 
' 12 d 24 d 
(3.9) 
which when inserted into the Bethe-Schwinger formula yields 
o(Ko—* zt 1 g , 
2(-! aT = Fir ).Uir- Fel¥o), (3.10) 
where ; 
yy ’ aye IM fa\* , 1 ,(a\*]-? 
Ui = U' NI +1 4 13 (ka)(4) nx W' —— oy 1-+nn).@'(5) | 
(3.11) 
We now restrict ourselves to the case when the ferrite sphere is situated 
at a magnetic antinode and the DC magnetic field is perpendicular to the 
rf magnetic field (Fig. 1). The element of Wire which contributes to the 
evaluation of (kj—k) is the (2,2) component of 
vw’ [T+ 4y’+o(I-+nn).p’]-, (3.12) 


3 
where _— alg) (3.13) 


x —j« 0 
and a = (i x 0). (3.14) 
0 0 0 


Assuming the approximate form 
= 4nM as 
Xe = AEA —je, = *" 
where for typical ferrites near the resonance DC field 
xX = «x = 4y.., we obtain to first order in a 


; 1+(4a)x. ae 
° = 4 .| —-— e 3.16 
(Better 3X Fesmsesis ( ) 
The imaginary part of this quantity has its maximum at a value of H 
given approximately by 
3 
at ) (3.17) 


A 
Huy = Hy + “2H 4 T(E 
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Fic. 1. Ferrite sphere near a wall. 
from which we find that 


Al = Mh m 0) Mil othe Se 


(3.18) 

Miles (6) has measured the wall effect for several ferrite spheres of 
different diameter and he obtained an average value of %M for AH,,.. 
Considering the effect of possible inhomogeneities in the sample, which 
would affect the experimentally determined value of M and also change 
the form of the elements of the magnetic permeability tensor, the agreement 
between our approximate theory and experiment may be considered fair. 
B. Size effect 

The quasi-stationary approximation which we considered in the last 
section shows no dependence on the size of the sphere. To obtain a size 
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dependent effect we must include higher orders in (ka). We therefore 

introduce the next approximation to the fields into variational expression 

(2.41). We set 

F(r) = A+B.r, F*(r) = At+B?.r, 

im 

where B.r | ; (3.19) 
Bur 

B¢ and B” being constant dyadics. 

Varying variational expression (2.41) respectively with respect to A+ 
and B+ yields the equations 


[ Adr = [ Ale)dr+k [[ Gir’ r).W'(A+B.4) drdr’ (3.20) 
and 
[ rB.rdr = [rA(r)drt+k ([ r'Gfr''r).W@'.[A+B.r] drdr’. (3.21) 


From the Bethe-Schwinger formula we note that for results for (ky—k)/k 
valid to the order (ka)3, we require A to order (ka)? and B to order (ka)?. 
But from equation (3.21) we find that B is at least of order (ka) and 
higher. We therefore approximate equations (3.20) and (3.21) by 


A = [1—1(k,a)?|Fo(t,) 4 am | | G(r'\r).W'.[A+B.r]drdr’ (3.22 


( rB.rd- | rF,(r) dr+k tf Grr). W@'.Adrdr’. (3.23) 
( rrdr = }(AV )a’I, (3.24) 
| | r’G,(r’ r) drdr’ (AV )o, F x Ta?n,(ka)j,(ka)], (3.25) 


and {fs S(v'or).D.rdrdr’ = (AV )o,{a*n,(ka)j,(ka)|D * F, (3.26) 
where D is a constant dyadic, and for the dyadics D, = ab, D, = cd, 
we define the operations 

D, «D, = (b.c)axd, 

D,:D, = (b.c)(a.d), 

D,.D, = (b.c)ad, 
and D, «D, = a(bx« c)d. 
We thus find that 


A = [1—Mkya)?|Fo(ty)— 4M’. A+ A( ka)?’ A <i M' xB (3.28) 
” 


and B! = V.F,(r)— ko, 1x (@’.A). (3.29) 
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Inserting equation (3.29) into equation (3.28) yields 


A= o>. 1—so(ho@)? |For) —" 0, W" X{VF (rp n'} (3.30) 


— 


where DG = F+4U'—s(kaPa’ +—— 0, (U"— FI U\)o,.U', (3.31) 


from which we obtain 
Bi — VF, (r,)—jko,F x (U' .G).F,(r,). (3.32) 


Wis the scalar of the dyadic matrix @’, which for wp of the form of 
equation (3.14) and € = ef is given by 
, Be’ 0 
Ww, = ("* (3.33) 
. 0 2x 
We now restrict ourselves to the case of a ferrite sphere in a rectangular 
cavity excited in a T'E,,,,, mode, which has only one component of electric 
field lying along the z-axis. For this case, since V.e, = 0, we find that 
Veg = &y X(T X Cp )€y €5 = jko(My  Cp)€o 6: (3.34) 
We are specifically interested in the case when the sphere is situated at 
either an electric or magnetic node as given by the equations 
sink,« = 1, and sink,y = 0or 1, 
where €y; = sink, xsink,y 
For the first case (sink, y = 9), 
“ = 0, (3.35) 
while for the latter Vh, = i #F ~k? ji+- k? ij). (3.36) 
0 
Performing the integrations in the Bethe-Schwinger formula (2.13), we 
find that 
2(ky—k) 
k AV 


For either case 


2 
— {1 ip(Ky @)* |F o( (To). WN + [VF 5( (To )]: |M’. B}. (3.37) 


A=Z-, li K(kya)? + Oe aX | F (To), (3.38) 


where X = (\ hy ) (3.39) 
0 € /€& 
while 
(UVF }):(@'.B) = Fol(ry).[_F#K+fkky XW .G").F (ro), (3.40) 
0 if e@,(ry) = 0 


where = 3X 2h? ke? , (3.41) 
k2 4x if hy(r9) - 
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Thus, to the order (ka)’, 


2(ko- k) at M , a 3.4? 
Ek AV = F(T). Ml .1¢- F (8 o), (3.42) 

, , 9 2k, ka* a ' ” @ ‘ « 
where Ve = V’ .B-* 1—$(kqa*) + = X}+13F Ka’. (3.43) 


For the case when the sphere is located near an electric node we seek 
the value of H,,, for which the imaginary part of (pir;),, i8 a Maximum, 


where 


, 12(ka)?_, . (ka)? ee ie 
toe = w[T+ gw — yy (20) |x 
45 45€, 


a) e€’ 


D]- J2 «’ * 
x [ : Higa)? eee | 4 Ke 
»” €o 


15 €o 


Near resonance we find that 


—1(kya)?+ 2k, kare’ | 
(the t)ee = Xs 1—3( Ee, (3.45) 
< + PX4 


where B = 4- a(ka)?*(6 + =) (3.46) 


€o 


and y. is given by equation (3.15). Im(pig,),, has its maximum at a value 
of H, given by 


0. 


(Hy)eos = (H- : 


45 


a) PS “(6 5 “\anM, (3.47) 


€o 


This agrees with the result obtained by Hurd (1), whose H, is our 
(H —47M), and whose H is our (Hy) ses 

When the sphere is near the magnetic node, the maximum absorption 
occurs at the values of H, for which im(€‘,,).. is a Maximum, where 


?/J-7)2\ @’ agp \2 fog 't -1 
Ett =e (; — -: So 2x.1)-£| " 


€o 3 45 Je& 45 Ho €o 


2h, ka? 


ab -4 kya)? + 2Hoh x] +408 Le)Xs. (3.48) 


0 
The imaginary part of (€,,,)., has a weak maximum at 

(Hi,).. = Hf. (3.49) 
We thus find that in moving the sphere from the electric node to the 
magnetic node, the resonance field increases by approximately 47M /3. 


This increase in the resonance field as we move from a region of zero 
electric field towards a magnetic node has been observed by Miles (6). 
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SUMMARY 
The solution of the linearized magnetohydrodynamie equations of a viscous fluid 
s generated in terms of one component electric and velocity vector potentials. 
\pplications of the approach to plane wave propagation and to a particular eigen- 
value problem are considered, 


1. Introduction 

IN a recent paper Ludford (1) has considered certain aspects of the problem 
of propagation of small magnetohydrodynamic disturbances in an inviscid 
fluid. In particular he has shown that the linearized equations split up 
into sets of four and five equations and that each set leads, on elimination 
of all but any one of the independent variables, to a partial differential 
equation of the same order. In the present note the separation of the 
linearized equations of motion of a viscous fluid is examined from the 
viewpoint of generating solutions in terms of one component vector 
potentials. It is shown that a general disturbance may be resolved into 
two parts. one generated by a velocity vector potential and the other by 
an electric vector potential. In both cases the vector potentials have only 
one component, that parallel to the permanent magnetic field, and the 
corresponding scalar functions satisfy fourth and fifth order differential 
equations respectively. For an inviscid fluid these reduce to equations 
obtained by Ludford. 


For plane waves the solution generated by the velocity vector potential 
becomes the v-mode solution obtained by Banos (2), and the other 
solution reduces to Banos’s p-mode solution. For the particular case of 


p-modes the present method provides a far simpler method of generating 
solutions than that of Banos. The present approach may also be used 
to study the propagation of magnetoelastic plane waves (1). 

The solution for standing waves in a fluid bounded by a rectangular 
box of perfectly conducting material is also obtained in terms of the above 
potentials. This problem has been solved by Ludford and it is shown that 
the present approach provides the solution more directly and in a more 
compact form. 


(Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 3, 1960] 
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2. General equations 
The linearized equations of magnetohydrodynamics are (1) 


J = ourlli—«?~, 
et 


divH = 0, 


cH 
—— = —curlE, 
Pat 


'y = E+pv~x Hy, 
Co 


: ; ae H, — § grad p+ ~ grad(div v)+vV2v, 
ct Po Po 3 


x + po divv = 0, (8) 


p = Ap. (7) 
H, E, v denote the deviations of the magnetic field, electric field, and 
fluid velocity from the unperturbed values (Hy, 0,0). The deviations of 
the pressure and density from the uniform state (po, p,) are denoted by 
p, p, and a, denotes the velocity of sound measured at py, pg. The macro- 
scopic constants » (magnetic permeability), « (dielectric susceptibility), 
o (electrical conductivity), v (coefficient of viscosity) are all assumed 
constant. 

By vector manipulation of the above equations it may be shown that 
the scalars p and divE satisfy independent homogeneous differential 
equations. It thus seems reasonable to seek solutions of the equations 
with either p or div E zero and hence divide the general solution into two 
classes generated by assuming either of these conditions. We consider 
first waves with zero excess pressure; such waves have been examined for 
plane-wave propagation by Banos (2) and, following his nomenclature, 
will be called v-modes. 

If, for a compressible fluid, p = 0, then (6) and (7) show that divv = 0 
and hence the v-modes represent possible solutions of the equations of 
motion for both compressible and incompressible fluids. 

If H is eliminated from (1), (3), and (4) then it may be shown that 


ae e 
a ee 
E at ;| “Pop 


Be) 
\|3 = — 1] ¥en S| Hy +m grad div(v x H,)+ 


ion . graddivE. (8) 
oa ct 


5092.51 
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From (1) and (4), 
—div p(v x H,) = f “x : | div E, 


0 


and hence, from (5), (8), and (9), we obtain 


R ; (V2 —e sal Lt = Po (- vV2}Vv + . grad p| 
| ét o ct? } o ct} p \\et Po | 


ec 


nfs? a 


|v < H,) x Hy) —H, x grad div(v x H,). (10) 


act 


If p = 0 then (5) shows that there may exist solutions with v.H, = 0 
and we may thus consider the class of solutions generated by v = curl(#H,). 
It is easily shown that 


| c ( l Ce ( >\p = oor Yaa | 
Be € | , V2 Po oe H? v2 n 
male unal)( iF sii ar i? 


“8 =0, (11) 


y" 


wis? : 


where the y-axis is taken to be parallel to Hy. For v = 0 equation (11) 
reduces to one obtained by Ludford from somewhat different con- 
siderations. 

Though a class of solutions of the equations of motion is generated by 
assuming the existence of the function ¢ it is clear that ¢ is not a 
generating function in the sense that all the components are expressible 
in terms of differential operators acting on it. The original equations 
show that the expressions for the field components in terms of ¢ are fairly 
complicated. If the conductivity is infinite, however, 


OM _ 7? cur(gli,). (1: 
i 


c cy 


E = wH, < curl(dH,), 


The forms of the field components in terms of ¢ are fairly simple for 
plane waves. In this case equation (11) reduces to a quadratic in the 
square of the wave number, and hence there are two possible plane wave 
solutions. This quadratic, with suitable change of notation, is identical 
with that obtained by Banos (3) for magnetoelastic plane waves. 

We now consider the class of solutions generated by div E = 0, and 
call this class of solutions pressure modes. The nomenclature is not im- 
mediately obvious but it will be shown later that the plane wave pressure 
modes obtained by Banos are generated by assuming div E = 0. There- 
fore, to conform with Banos, it seems reasonable to call this class of 
solution pressure modes. 
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3. Solutions generated by a solenoidal electric field 
From (1), (3), (4) it is seen that 
| se a. i - J.H, = _ hy fy div E. 
| op ot a ct pe cy 
Equation (13) thus shows the possible existence of solutions with 
E.H, = J.H, = 0 
and we consider the class of solutions generated by E = curl(Hy,). 
We consider the general case of a compressible fluid, the corresponding 
results for an incompressible one may be deduced by letting a) > x. 
From equations (5) to (7) 
nite mm ae x H, +43 grad(div v)+- ’ © grad div v t wee, (14) 
ot? Po et 3 ct ct 


and hence 


Pols (+9 | laivy - 7 = |v {p= «| Ho. curl E. 
¢ 


Then, from (14) and (15), 


je Bilas 7 a oe | 
Polar “OT 3S et , (ae at\” 


ee 2, Wve yld H+ 
Pica \OT 3 al ja * 


brad (ai +5 +4 + ‘) |v an ‘a|Ho-curl E. (16) 


Equation (4) now gives 


(Pao) -S mm 


pe ee 4 é - 
= —*—H, x grad c i+ = — || ran =| H,.curlE. (17) 
Po 3 et|| ct 


In terms of the scalar function % defined by E = curl(H,), (17) becomes 


aa (8+ Falla ead) a) a = ea) 


BAG 2 2 S| 
oo (tall zal [P ale 


For v = 0 equation (18) was obtained by Ludford from entirely different 
considerations. For plane waves in a viscous medium there will be three 
possible wave numbers and, with suitable transformations, these may be 
obtained from Banos’s results for magnetoelastic waves. 
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For general time behaviour it is clear that % is not particularly useful 
for generating the field components. If we assume a time factor exp iwt, 
however, then 

;, l ' 
piwH curl(curl {H,), pv, = {_ [V2+ uw? |] — 1\ grad, x, 
\iwpo j 
(19) 
where the suffix t indicates the component perpendicular to Hy. For an 
incompressible fluid v, is obtained from div v = 0, and for a compressible 
fluid it is given by 


\ 


} : iawV2ln, 


{ lw 
2 2 
HPo t (a3 T3 


cy” 


oy ” 
(ag = \( _we4 cust) —1)(¥* S)e. 
3 /\twyo cy} cy 


For an incompressible fluid equation (18) becomes 
é a” Ree e2\ eV? Hi 2 
pa ~ Lede 2) p82 

op ot* ct Po Cy? 

The original system of nine first-order equations in time has now been 
reduced to the fourth-order equation (11) and the fifth-order equation (18). 
Thus a general solution may be effected in terms of the functions ¢ and y. 


\ 
ct 


hee (21) 


These functions are particularly useful for investigating time harmonic 
cylindrical waves and plane waves. For the particular case of plane waves 
the use of % provides a far simpler and clearer method of generating the 
pressure waves of a compressible fluid than that of Banos. This latter 
method consists of assuming a form for the velocity vector in terms of 
an angular parameter y which is determined once the wave number has 
been obtained. The use of % avoids introducing the unnecessary para- 
meter y and the solution may be obtained directly from y. 

We now apply the present approach to a particular boundary-value 
problem considered by Ludford and, following him, assume that 


yon o-! = ¢ = ©, 


4. Standing waves in a rectangular box 
The gas is assumed to fill a rectangular cavity (x = 0,a; y = 0,6; 
z = 0,c) in a perfectly conducting rigid body into which is frozen a uni- 
form magnetic field, the corresponding external field in the cavity is taken 
to be H, along the y-axis. The boundary conditions are 
=0 onz= 0,a, 
=0 ony = 0,b, 


»-=90 onz=29,C. 
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For v-modes we have 


¢ = | Fil 2)eos St +. F(x, z)sin - “t [expC 
Ay Ay 
where A, is the Alfvén velocity and F,, F, are arbitrary functions. From 
(12) it is seen that all the conditions on y = 0, 6 are satisfied if 4 vanishes 
on these surfaces; hence F, = 0 and w/Ay = mz/b, where m is an integer. 
The conditions on the other surfaces are seen to hold if 4, and hence F,, 
vanishes on these surfaces. 
For pressure modes equation (19) shows that all the conditions on the 
planes parallel to the y-axis are satisfied by % = f(y)cos(/z2/a)cos(nzz/c), 
where /, n are integers. From (20) it is seen that f satisfies 


fo4 BE’ + Cf = 
where Ba (AG+a5) | (7 ‘+ 4 


a2 Az ee ¢ 


and AzazC = wi—(ai+-A§)w?2*(l?/a? +-n?/c?). 


The electric field components vanish on y = 0, 6 if f(y) vanishes on these 
surfaces. It can be shown from (20) and (22) that v, = 0 on y = 0, b if 


Gf ip ff 
dy’ = ga (a0 B wT, 


- ony = 0,b. 
dy 


If it is now assumed that f is proportional to exp(iky), then a transcen- 
dental equation is obtained for k, which was first obtained by Ludford; 
his approach is rather more involved and requires the introduction of an 
additional function g(y) and of differential equations relating f and g. 
It can be shown that these equations reduce to (22) and that Ludford’s 
boundary conditions are identical with the above. 
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SUMMARY 


excitation of surface waves by a plane electromagnetic wave incident on a 

dye is investigated. It is assumed that the impedances of the two wedge surfaces 
onstant and the mathematical problem is the solution of Helmholtz’s equation 
utside a wedge on whose surfaces impedance conditions are satisfied. If one 
impedance is zero then it 1s shown that, provided the other impedance is such that 
it supports surface waves, these will be excited unless the ratio (exterior wedge 


ingle) wis | 2n, where n is an integer. If the impedances of both surfaces are non- 


zero and independent then surface waves are excited except when the above ratio 
is 1/n. For only one impedance non-zero explicit forms are deduced for the amplitude 
of the surface wave for large and small values of the impedance. The solution may 
also be employed for the diffraction of an H-polarized plane wave by a metallic 


wedge 


1. Introduction 

[vy the impedance of a plane is such that it supports the propagation of 
surface waves then it is known that these waves will be excited when 
certain electromagnetic fields are incident on the plane. A plane wave 
does not produce a surface wave whereas the field of a line source does. 
In view of the fact that sharp edges, when excited by an arbitrary field, 
behave as line sources it seems likely that a surface wave can be produced 
by a plane wave incident on a wedge whose surfaces support the propaga- 
tion of surface waves. It is thus desirable to investigate the conditions 
under which surface waves may be excited by a plane wave incident on a 
wedge of arbitrary angle. This is one object of the present investigation. 
In order to achieve this it is necessary to solve the problem of the diffraction 
of a plane wave by a wedge of arbitrary angle with different, but constant, 
impedances prescribed on its two surfaces. A solution of the diffraction 
problem for a right-angled wedge with zero impedance on one surface is 
given in (1). 

The approach of (1) cannot be generalized to cover the case of an arbitrary 
wedge and the method developed by the author for diffraction by an 
imperfectly conducting wedge (2) has to be employed. One particular 
result that may be deduced from the general solution presented here is the 


(Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 3, 1969) 
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solution for the diffraction of an H-polarized plane wave by an imperfectly 
conducting wedge. This solution may be deduced by suitable changes of 
notation from (2) but its form is such that considerable manipulation is 
necessary to demonstrate its behaviour in the limiting case of infinite 
conductivity. The problem of diffraction by an imperfectly conducting 
wedge has also been solved independently by Senior (3). 

It is shown that, if the surfaces are such that they support surface 
waves, these waves will be generated under most circumstances. If one 
surface impedance is zero then the surface waves will not be excited if the 
exterior wedge angle is an integral submultiple of $7. For equal surface 
impedances surface waves are not excited if the exterior wedge angle is an 
integral submultiple of 7. It is also shown that there may exist a relation 
between the surface impedances which is such that no surface wave will be 
propagated along one of the surfaces. 

For the particular case of zero impedance on one surface expressions are 
obtained, for large and small values of the impedance, for the amplitude 
of the wave propagated along the other surface. For large values of the 
impedance and a right-angled wedge agreement is obtained with (1), for 
small impedance the results do not agree. The reason for this discrepancy 
is that, for small impedance, the expressions in (1) are incomplete and 
terms have been incorrectly neglected. 


2. Statement of the problem 


A cylindrical polar coordinate system (p, @, z) is chosen with its axis along 


the vertex of the wedge and such that the wedge occupies the region 


y<O0< Qn. 
The general mathematical problem to be considered is the solution of 


(V?+k*)p = 0, (1) 
subject to the conditions 


cp 
a 


ikn, on 6 0, % —itkn, on 6 Y; 

co > 
where & is real and », and », are constants assumed to be in the first 
quadrant of the complex plane. It is also assumed that a plane wave 
exp tkpcos(@—6,) is incident on the wedge. Thus it is necessary that 4, 
apart from the terms of geometrical optics, satisfies a radiation condition. 
We shall employ the radiation condition lim p*(é¢/ép+ikdé) = 0. In 


| 
terms of a solution of the time harmonic wave equation this corresponds 
to assuming a time behaviour exp iwt, where w is positive. 
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The solution to the boundary value problem may be deduced by a 
generalization of the results of (2) and we have 


°. 


-f [aie +o) ee It Gry M—ik ; 
d || (w+-@) io 5-=. (w—O@)|exp(—ikpcosw) dw, (2) 


c 
where 
G(w) G(w T 2y)K(w, 1 )A(w, ies 1, )K(w+y, a) K(w+y, i Og), (3) 
and K(w,«) = tan }(w—42—a), 
where «,, a, are defined by cosa, = »,. Our assumptions about 7, show 
that re x, < 47. The contour C is taken to be below all the singularities 
of the integrand and with its end-points at —io0oo—4m and —ioo+ $n. The 
integrand is required to be bounded at the end-points of C and to possess 
no poles in the strip im w < 0,0 < rew < =. 
If H(w) denotes a particular solution of (3) then the general solution is 
G(w) H(w)H,(w), where 
H,(w) = H,(w+-2y). (4) 
By arguments employed in (2) it can be shown that the conditions at 
infinity are satisfied if the integrand of (2) has no poles in the lower half- 
plane with real parts between 0 and z and if 
G(w)|sin w+ y,| = —G(227—w)[sin w— y, |. (5) 
If », and », are both set equal to zero then (3) and (4) are identical and 
5) be 48 . . 
(9) become G(w) = —G(227—w). (6) 
Thus, if ,(w) is chosen to satisfy (4) and (6), H(w) may be chosen so that it 
tends to unity as the impedance tends to zero and 


[sin w+», |H(w) = H(27—w)[sin w—y,]. (7) 


This method of splitting up the solution is particularly suitable for small 
values of the impedance. Equations similar to (6) and (7) were derived 
in (2), in this case the sign of H on the right-hand sides was reversed and the 


solution thus obtained is more suitable for large values of the impedance. 
Also, from (4) and (6), 


Hj(a+ry)=90 (r= 0, +1,..., +00). (8) 


3. Solution of the difference equation 


A solution H(w) may be obtained in terms of the function F(w) obtained 
in (2). We shall summarize the necessary properties of F(w). It is an 
analytic function of w apart from poles on the real axis at 


w= —2n07—2ry— tr 


and at w = 474+2(r+1)y+2nz, n, r = 0, 1,..., 00. It also has zeros at 
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w = 2(r+1)y+2(n+1)r—4n and at w = —2(n+1)r7+4n—2ry, where 
n,r = 0, 1, 2,..., 00. For large values of w 

F(w) ~ (27) exp}i(rw—7) = (imw > 0), (9) 
F(w) ~ (27) exp{—}i(rw—7)} (im w < 0), (10) 

where + = 7/y. F(w) also satisfies the equations 
F(w) = F(2y—vw), (11) 
F(w)F(w—) = 2i cosec 47(w—4n), (12) 
F(w) = —tan }(w—42)F(w+2y). (13) 
Clearly a solution of (3) which has the appropriate behaviour for zero 


impedance is 
F(w—a)F(wt+y—ay) 
F(w—7+-0,)F(w+y—7+ a.) 





H,(w) = (14) 


From (11) and (12) we obtain, after some manipulation, 


Hw) _ [sinw—n,]|sin $7(w+a,— $2) cos $7(w+ a,— $7) 
H,(27—w) — [sin w+-n,]sin 47(w—42—a,) cos $7(w—47—ay) 





(15) 


Thus, from (15), a possible solution of (3) and (7) is 
H(w) = 4H,(w) sin $7(w— $2—.«,) cos $7(w— 42—a) cosec r(w—7), (16) 


and H(w) tends to unity as the impedances tend to zero. 

It now remains to determine H,(w). This function has to be chosen so 
that the integrand of (2) satisfies all the conditions imposed on it. H(w-+6@) 
is bounded at the points of C and, hence, so must H,(w-+-@). The integrand 
is required to have no poles in the strip im w < 0, 0 < rew < 72, apart 
from those which produce the geometrical-optics terms. H(w+-@) may 
possess inadmissible poles and H,(w-+6) may thus be required to have zeros 
which cancel these poles. The properties of F(w) show that the only 
inadmissible poles are at w = 7+6-+ry and, from (8), H,(w--@) will auto- 
matically have zeros which cancel these poles. H(w) is independent of 6, 
and hence the geometrical optics terms must be produced by H,(w+-@). 
It is easily verified that the solution of (4) and (6) which produces the given 
incident wave is 


H,(w) = [cot 47(w—a—6,)+cot $7(w+0,—7)]. (17) 


l 
2iyH (+65) 


4. Behaviour at infinity 

The behaviour of ¢ at infinity is obtained by deforming C into the 
saddle-point curves through w = 0 and 7. These curves are illustrated 
in (2). Some poles of the integrand may be captured in this deformation 
and this part of ¢ will be denoted by ¢,,,. The integral over the saddle- 
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point curves is essentially a diffracted field and will be denoted by ¢dyiq. 
By a transformation similar to that of (2) it may be verified that 


duit | | H(w+-0)H,(w+6)— 
S(0) 
H(w+-0+ 27) H,(w+-6+ 27) lexp(—ikpcosw) dw, (18) 
where S(0) is the saddle-point curve through w = 0. 

For most values of 6, dai may be evaluated by the normal saddle- 
point method. If there are poles of the integrand near the saddle point, 
however, the modified saddle-point method of Clemmow (4) and van der 
Waerden (5) must be employed. The method is straightforward, the 
calculation and the resulting complicated expressions will not be repro- 
duced here. 

Our main interest in the present analysis is in the terms obtained by 
capturing poles of the integrand of (2). It is easily verified that the poles 
of H/,(w) produce the plane wave terms of geometrical optics. It therefore 
remains for us to evaluate the contribution of the poles of H(w) to 4,,.. 

The properties of F'(w) show that there are poles of the integrand of (2) at 
w kara, —0—2ryandatw har +-ag+(2r—1)y+0,r = 0,1, 2,..., 00. 
We consider the contribution of the first class of poles with r — 0; it is seen 
that it gives rise to a term which behaves as exp|ik{ny—av(1— ?)}], which 
is of the surface wave type. Thus it is seen that there is a possibility of 
obtaining surface wave terms in the present problem. The second class of 
poles clearly produce surface waves on the surface 6 = y. For simplicity 
we consider only the first class of poles ; this corresponds to investigating 
the surface waves on 6 — 0. The analysis for the other surface may be 
carried out similarly. The analysis is further restricted to the case r = 0; 
it may be seen that the higher values of r correspond to multiply reflected 
surface waves and the properties of these may be obtained by a simple 
extension of the following analysis. 

The conditions under which the pole producing the surface wave term is 
captured will now be obtained. Since the real and imaginary parts of 
7, are both positive it is seen that «a, = a,—tb,, where a,, b, > 0. Itis now 
easily verified that the surface wave term is captured for 


6 < a,—4n-+-cos~\(sech b,). 


This is also the condition for the propagation of a surface wave along a 
plane surface when a line source is placed at the origin (6). Thus it seems 
that the vertex of the wedge effectively acts as a line source as far as the 
production of the surface wave terms is concerned. If 


42—a, > cos~!(sech b,) 
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then no surface wave will be produced. If the impedance is of the type 
produced by a highly conducting surface it may be shown that the surface 
wave is not produced (cf. Clemmow 7). For small values of », the pole of 
H(w +6) will be near the saddle point and the modified saddle-point method 
will also have to be employed in this case. 

It may be shown that the amplitude of the surface wave is 


[—7i sin $27 cos $7(a, —ag—7) F(2y— $a) x 
< Fly+ %y — Me ~ ha) F(x + Og +y+ x a) Hol x 


C (19) 


7)| 
sin r(a, $m) F(y+a, +- My — 3m) F(7-+-65+ y— a) F(2a,- a 


We now examine the conditions under which C, vanishes. The zeros of 
F(w) are clearly all on the real axis and thus, apart from the zeros of 
sin $a7, zeros will only occur if ima, = 0 or im(a,—a,) = 0. The above 
analysis, together with the fact that a, < $2, shows that the surface wave is 
not captured if ima, — 0; thus we require im(x,—.,) = 0. Unless a, = a, 
this relation will only hold if a fairly complicated relationship holds between 
7, and »,. We shall not consider this special case in any more detail but 
merely point out that it may be possible to choose 7,, 7, so that the surface 
wave on @ = 0 is not excited. Thus, for unequal impedances, C, — 0 if 
y 7 2n and for equal impedances C, vanishes for y — 7”, n being an 
integer. If », = 0 then 
een 2a sin $7 Hy(a,— $7) F(2y ~hm) (20) 
sin 47(x, — $7) F (2x, — 37) 
The above expression is somewhat simpler than that of equation (19) 
but it is still formidable and it is thus of interest to examine whether any 
further reduction is possible. The two cases of interest are large and small 
values of the impedance, these correspond to values of a, near —ix and 
sm respectively. The case of values of a, near —i%x is considered first, 
After some manipulation it may be shown from equations (9), (10), and 
(20) that 
0 


A exp({irz)sin $77 sin $76,(27,)-*7 F(2y—4z). (21) 


Sr 
~ Qa )' 
For the particular case of y = 3a considered in (1), it has been shown 
(2) that F(w) = (27)! cos }wsec }(w—4z), 
and thus, in this case, 

C, = 4(2n,)-*sin 44, exp(jiz). 
The value of C, for large 7, has also been considered in (1). In comparing 
the present results with those of (1) it should be remembered that the 
solution in (1) satisfies the radiation condition lim pt(ef ép—ikf) — 0 and 


p-> 
it is assumed that », , m. ~ 9. If these dethrone and the notational 
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differences are taken into account then it may show that the results are in 
agreement, 

For small », we have «, = 47—7, and thus, from (13), it may be shown 


that C, & 271, H,(0) 


77, 8iN r77/COS 7A, — COS Tz. 


This result does not agree with the corresponding result in (1). The dis- 
crepancy is due to the fact that, for small 7,, the expansion obtained in 
(1) for C, is incomplete. In the analysis of (1) several expressions of the 
form {(7,)+-4,g(m,) occur and in all cases these have been calculated by 


replacing f, g by the first term in the Taylor expansion and non-zero terms 
of the form »,f‘(0) have been neglected. The analysis of (1) is correct only 
to the zeroth order and it is shown that the zeroth order term in C, vanishes. 

It should be noted that the infinity in (22) for 4, 7 is due to the fact 
that the expression for C, has been obtained by expanding H,(7,) about 
», — % For 6) near 7 this expansion is invalid as H,(0) becomes infinite 
for 6, — 7. Thus, for 6, = 7, H,(0) in (22) should be replaced by H,(— )). 
The resulting expression will not be infinite for 6, = 7 as H,(—7,) only has 
poles for real values of «,. Such values of «,, however, are excluded as they 
do not produce surface waves. 

The solution from an H-polarized wave is obtained by putting 
ny = Ne = 7H. Where 7 is defined in (2). 

The form of 1/(w) for small 7 shows that it may be expressed as a power 
series in 7 and a solution could be obtained by a perturbation approach. 
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SUMMARY 

The equation governing the hydromagnetic oscillations of a liquid globe rotating 
in a uniform magnetic field is solved in cylindrical polar coordinates and the solution 
adapted to satisfy spherical boundary conditions, The period equation is obtained 
in the form of an infinite determinant. Numerical approximations in the case when 
the Coriolis force appreciably affects the motion indicate that convergence is likely 
and confirm the order of magnitude of the period as obtained by Cowling (1) for 
the case of a liquid and a gaseous globe. 


1. Introduction 
Tue problem of toroidal oscillations of a rotating fluid sphere is of interest 


in connexion with theories of the solar cycle of twenty-two years and with 
the problem of the secular variation of the earth’s magnetic field. 

Torsional oscillations were first discussed by Walén (2) with reference 
to the solar cycle, while Plumpton and Ferraro (3) found the first exact 
solution for torsional oscillations in non-rotating stars and showed that 
these were such that shells generated by rotating a magnetic line of force 
round the axis rotated independently. Moreover, because there is a con 
tinuous set of such shells there is also a continuous set of free periods, an 
unusual result. Cowling suggested that a discrete set of periods would be 
obtained if the fluid globe were rotating and showed that the introduction 
of the Coriolis force would slow down the oscillations considerably in the 
case when the globe was a uniform liquid. If the globe were gaseous he 
showed that the effect of the Coriolis terms would be small and conjectured 
that the effect of the rotation would be to select a discrete set of periods from 
the continuous set of Plumpton and Ferraro. 

In this paper we have attempted to give a more complete discussion of 
the oscillations of a rotating liquid globe. The equations of the problem 
were first derived by Cowling (1) but not solved by him; they were later 
derived also by Namikawa (4) who considered the solution of the equations 
for a liquid cylinder rotating about its axis. Namikawa showed that for 
stationary waves of given wavelength parallel to its axis, the cylinder 
has two natural periods of oscillation. As the Coriolis force becomes large 
the longer period corresponds to the hydromagnetic modes of oscillation and 


(Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 3, 1960) 
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the shorter period to the mechanical modes. He finds the period of hydro- 
magnetic oscillations to be 7-7 x 10° years for a cylinder whose dimensions 
are of the same order of magnitude as those of the earth’s core. 

The solution given in this’ paper will be restricted to the case of axial 
symmetry, although the non-symmetrical oscillations may be similarly 
treated. Since spherical boundary conditions are imposed upon a solution 
which is essentially cylindrical in character the problem may be of some 
general interest. 

Our results for the axially symmetric case are in agreement with the order 
of magnitude derived by Cowling for the periods of free oscillation for the 
fundamental modes though we have only been able to obtain a fair approxi- 
mation to these periods. It is hoped to improve the approximation by the 
use of a digital computer. As Cowling pointed out, the results obtained 
cannot be applied to real stars; however, the core of the earth is generally 
assumed to be in the liquid state and the results are of some interest in 
connexion with possible existence of free oscillations of the fluid core in the 
permanent magnetic field of the earth. Assuming infinite conductivity we 
find that the free periods of these oscillations are of the order of tens of 
thousands of years and too long compared with changes of the earth’s 
field to be of interest for the secular variations. 

The conductivity may affect these oscillations in the earth’s core; in fact, 
as has been pointed out by Chandrasekhar (5), if one adopts the values 
given by Bullard and Gellman (6) for the conductivity, namely, 3 x 10-6 
e.m.u., the decay time for the earth’s field is of the order of 14,000 years. 
This would be comparable with the periods of free oscillations. But since the 
free periods of hydromagnetic oscillations themselves are about 100 times 
larger than the period of time of 100 years during which the earth’s field 
changes appreciably, this would suggest that the secular variation of the 
earth's field must be the result of adynamo mechanism due to motion in the 
core across the general magnetic field pervading it. 


2. The equations of the problem 

We consider a liquid sphere rotating with uniform angular velocity w, 
in a uniform magnetic field Hy parallel to the axis of rotation of the sphere; 
the angular velocity w, is supposed, moreover, to be so small that we may 
neglect the departure of the surface from the spherical form. We shall 
further assume that the conductivity of the liquid is infinite. Choose 
cylindrical polar coordinates (aw, 4, z) such that the z-axis coincides with the 
axis of rotation and is directed in the sense ofw,. Denote the density of the 
liquid and the pressure at any point by p and p respectively and by (u, v, w) 
the resolutes of the velocity v of the liquid in the directions of w, 4, z 
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increasing. Likewise let (h,, h., h,) be the resolutes of the small variation h 
in the magnetic field and denote by 2 the gravitational potential. Let the 
suffix 0 refer to steady state values and primes to small variations from 
these. The linearized equations of motion in cylindrical polar coordinates 


are then 
_ ee op’ eQ’ = Ay[eh, ch, 
Pol - — 2w wow — + —py—- + — 4}, 
ct Cw Cw A4n\cz Com 
ev’ H,, eh, 
, + 2u’ = —? ° 
( et us) 4x éz’ 
eu’ Cp’ eQ’ 


Po z= = coe (3) 
ct C2 


where we have written v’ = ww’. Since the conductivity is infinite the 
lines of force are frozen in the fluid; the linearized equation expressing this is 


curl(v ( H,), (4) 


the resolutes of which are 


ch, eu’ ch, ev’ ch,  —H, @ 


H,-—, —- = H,—, — = ——(wu’). (5) 
ot oz ot oz ct aw Cw 
The equation of continuity of the liquid divv = 0 and the equation 
divh = 0 can be satisfied by introducing the liquid and magnetic stream- 
functions y and L’ respectively, such that 
wen ®, (6) 
Cz Cw 
eu 


ou 
hw = - ’ h,w =; . 
Cz Cw 


uo = 


(7) 


Substituting these expressions for the velocity and magnetic field resolutes 
in (1) it becomes 


a czcl Ow 


eat ower = (Psa) 4% AU, (8) 


Po Amp) a 


where we have written 


ofl a oe 
A rod ( )+ ; 


éw\w ew) | ez? 
likewise substituting (6) and (7) in (2) and (3) we find 
Cw’ —- 2wo, Cab _ H, chs 
et ow Oz Ampy mw Cz 


and 1’ = - bd (2 +0") 


w Cac 02\ po 
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Eliminating ( p’/p»)+-2’ between (8) and (11) we obtain 
é ‘ a) n Pn H, (=) 
Clim C2 4np, w C2 
Again from (5) we have 
ch. ov’ Ca" 


Hy - oH, 


ct Cz 02 


(13) 


and substituting in (10) after differentiating with respect to time we have 
“2 a2 2 2, 
- Qt a (14) 
et? C2? uw” czct 
where }, is the Alfvén velocity defined by 
I H,/(47po)*. (15) 
The first and last of equations (5) give 
ay ' 
: —. (16) 
ct Cz 
and, using (12) this becomes 
2 ne 25° 
(' .—¥5- ,) AY + Levy 2 = 0. (17) 
ct* 02° czcl 
Eliminating ys between (14) and (17) we obtain an equation for w’ in the form 


o O4u,' 


vs Ss } Aw’ + deve or? 0. (18) 


et? é2? e2*ct® 

These equations were first derived by Cowling (1) and later by Namikawa. 
Gravitational forces do not appear explicitly in (18) since, as was pointed 
out by Cowling, a uniform star suffers no change in gravitational energy 
after an internal displacement. When w, = 0, equation (14) reduces to 


the equation a 
C°wW 


; (19) 
ct 
derived by Plumpton and Ferraro. The absence of the coordinate w from 
this equation together with the appropriate boundary conditions indicates 
that the sphere rotates in shells generated by rotating a line of magnetic 
force about the axis of symmetry. 


3. Solution of the differential equation 
We assume harmonic variations of period 27/c, so that w’ x e'; sub 
stituting in (18) this becomes 


{ e2 \2 e2 
| (o? T Ve 3) $i? o | , 
cz" C2 
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Let us seek a solution of this equation of the form 


; : = 
w = a Pa xX, 7) a), (21) 


where a,(z) are functions of z to be determined, a,, ay,..., are the zeros of 
Bessel functions J, of order 1, and a is the radius of the sphere. Substituting 
21) in (20) we find 


d2\2/d2 2 2 
(02 V2 a, (7 =) 42 o? , 


2) ‘ 
| 9 122 dz? a? a,,(z) 0, 


f 
| dz* 
a linear equation with constant coefficients, the solution of which is of the 


form ¥ A, exp(i8,,z/a), where the A, are constants and £, are the roots 


_— ” 


of the characteristic equation 
(a*a? — V3 B? )?(B2 + a2) — 4w2 ato*p? 0. (23) 
By considering the sign changes of the left-hand side of this cubic in p? 


2 where A,,, 4, v,, are real and 


, 
ne ‘ 


we find that the roots are B2 — A®, w2, —1 


n 


such that 
L - 2 =< —a*o?/V?2 - ge<O< ak << mw. (24) 
The solution of (20) is therefore of the form 


oe 
i ¥ {A,, cos(A, 2/4) + B, cos(p,, z/a)-+C,, eosh(v,, za) +4 
n” 1 


FE, sin(A,, 2a) + F, 


n " 


sin(y,, 2 a)+G, sinh(v, 2 a)}J,(4, wa). (25) 


This can be written in the abbreviated form 


wu’ 2, r b 2 [A,, cos(B,,z/a)+ EF, sin(B,,2 a)|J,(a, aa) (26) 
n Ba niftn Vn 


provided that we replace G@, by G, i inthe sum. The presence of the hyper- 
bolic terms in this solution is analogous to the skin-depth found by 
Roberts (7) for the reflection of Alfvén waves at a plane boundary separating 
two media. 


4. The boundary conditions 

At the free surface of the sphere the boundary conditions to be satisfied 
are that the excess pressure shall vanish there, that is p’ — 0 on r =a, 
and that the three resolutes of the magnetic field shall be continuous. The 
latter condition also ensures that the electric field resolutes are continuous 
at the boundary; for the electric field outside the sphere consists of the 
sum of a meridional part derived from a scalar potential and an azimuthal 
part proportional to the magnetic stream-function outside the sphere. 

02 51 v 
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Since the azimuthal electric and magnetic fields inside the sphere are con- 
nected by the same relations, it follows that the azimuthal electric field is 
continuous at the boundary. The meridional field outside the sphere is 
determined by the boundary values of the tangential resolutes. 

The method of solution consists of expanding the fields in terms of 
Legendre polynomials; in practice it is simpler to multiply the series by the 
appropriate polynomials and integrate over the boundary. In this con- 
nexion the following integrals are used and are quoted here for future 
reference. 

- | : oe 
elecostoos’ J (z sin @sin &)P,(cos @)sin 6 dé = (= 


}PiP,(cos P)Jp.4(z), (27) 


. ») 
elzoosFoosy J (z sin Osin ys) P;.(cos @)sin?0 dé = (=)! ik-1 P’.(cos wh), , ,(z)sin ys, 


0 = 

(28) 
where J, and J, are Bessel functions of order zero and unity, P, and P;, 
are the Legendre polynomial of degree k and its derivative respectively 
and J,,., is a spherical Bessel function. 


5. The pressure equation 
We must first consider the variations in the gravitational potential 
brought about by a distortion of the surface. Assume that this has the form 


r=a(l-4 2 (met), (29) 
mam 


where e, is small, 4 = cos @, and @ is the colatitude of a point on the sphere. 
Then we find that at the surface 


srGpa® SY {(k—1)/(2k+ lfe, P.(u), (30) 
ko1 
where G is the constant of gravitation. 
Again from (6) the radial displacement &, of a point on the surface is 
Ow "Cw C2 


(u' sin 6+ w’ cos @) ico = — - (" cos 6 “sin®). (31) 


From (14) and (26) it then follows that 


io — < aigtQ— V2 R2 Bz B,,z = 
i) . ~ ’ “ € 
% = > * o"| A, sin—"- — E,, cos )a( } (32) 
a 


2w9o — an ap, a a 
n Br 


omitting the time factor in this and all subsequent equations. Thus from 
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(31) we have 


l . ree 4 x ‘x, 7). 2 
1S (a2o?— V2 B2)| —" J,{ “— sin B, cos 6— 
» ae" OMn 0 
2w,o7a?| B,, a a 
n 1 Bn 
zx 


Jf" }ooa(s=)ein 44+ > > otot—VaRt)> 


yee *n sin Bn cos 6—J, Xn sin Bn sin @ E,). (33) 
B, “\ @ a a a / 
Multiplying this equation by P,(cos @)sin 6 and integrating from: @ = 0 to 
6 — mand using (27) and (28) we find after some reduction that 


79 99. 


2aek “<< [ a a’o?— Vi B (77) , 
— = - 2 3 ___ Se. Pi seond. Wi. .slev,) x 
2k+1 3 “ a 3.)( 2w, 07a” = ut $$) . - n) 
x[A, rei*+E,imi*]. (34) 


In this equation we have introduced an amplitude w, and a phase ¢,, 


suc é i 5 
such that B,, = w,, cos ¢,,, x, = w, sing, (35) 


and that when we set 8, equal to X,,, w,,, tv, in turn, so we let w, = 2,,, Yns Zp 
and ¢, = &,, ,, ¢, in turn, where 
z, conf. =A x, sing, = a,, (36) 


Y, COS, = Has Yn SIN Ny = My, (37) 
z, cost, = t,, z, sin ¢,, = ; (38) 
Thus z, = i(v?—a®)! and, since v, > a, by (24), it follows that cosZ,, is 
real, though not €,, itself. 
Multiplying (30) by #,sin@ and integrating from 6 = 0 to 7 we have, 
using (27) and (28), 
] b.. fi o x, (a2? — V2 82) 
Q P,sin@ dé = —8nGp|— S 7 En Ex 
| oe “7 | 2 — r 2w, ow, aB,, 
n ” 


x P(cos d, )J,.,(w, [A, retk¥+£, imi*] (k = 1,2,3,...). (39) 
From (11) and (14) we obtain 
a2 fy! 4 ( 2 [22 Ad 
: it ie i ~¥— 
€2°\ po ow Cm|2w,\ct? ez? 


whence, using (26) and integrating, we have 


a” 


(40) 


72 Q2 


i S > x, {(a2o?— V2 B2) A cos(P = my 
2 n a 


» << 
_ aw | 
n=1 Bp 0 n 


+E, sin(* *)|}o(™2 | arbitrary constant. 
a a 





292 P. C. KENDALL 


At the boundary p’ = 0, and using (27), after multiplication by 7, sin @ and 
integrating from @ = 0 to 7, we have for k = 1, 2, 3... 

Pts <. 5 x, (a20?— V2 B2)(27\4 

| Q’P. sin 6 dé S S : aaa a } 1A(COS p,, ey .4(W,) » 


/ » 2 
mpm - 2w api, 


0 n 1 j n 


<[A, rei*+-E,imi*l, (42) 


The case k = 0 merely serves to evaluate the arbitrary constant. From 
(39) and (42) we have then the infinite set of equations for the unknown 
constants A,, , ae 


n m seen 


ifs) (2 }(/P.(cosd,)  87G@p/k—1)\ P;(cos¢,,)) 
}" B,, 30° 4 u | 


w 


” = 


n 


(wv, [A, rei*+E, imi*] = 0 (k 


This infinite set of equations ensures that p’ = 0 on the boundary r 


6. The magnetic field variations 
We next consider the field variations; from (14) and (16) we find 


ae Vic a)e (44) 


» 2 2 
~-Wy oO 


On using (26) we have for the magnetic stream-function inside the sphere 


ol, SS (ao? v363) A. cos( ‘} +E, sin(* ‘| -( =n ”), 
a a t 


» 27,2 
=-W, 7 C1 hme oe ( 
n 


(45) 
from which we derive the variations in the field resolutes as 


My ~ » B,,(a?o? v3) A © sin(n*) 
, a 


2w, 0a —— 
n l 


Oe cos( = 
SS etter ViR%)| 4, cos(” 


2w, o7a® a } 
4 i li x, W 
| E,, sin hn Je —n 
a a 


1 
Outside the sphere (neglecting displacement currents) the magnetic 
stream-function is of the form 


Br 





U’ p D,, r-"a"*+! Pl (m),/(1—p?), (47) 


n 1 
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whence the field resolutes are 


a ae ->3 n a"*In Ph, (1), (48) 


a OZ 


» 180 — D , 
a — n ne l 
hy oe 2, pit a" *In(n + L)P, ..(H). (49) 
Continuity of the magnetic field at the surface of the sphere implies that 
h, = h, and h, = hj; multiplying these two equations by (1—*)!P)(y) 
and (1—y?*)!P.(u) respectively and integrating from @ = 0 to 7, using 
(27) and (28), we find 


= aH, = 2a Rv, 
Sw o%ad ae a, (a2o? — V2 B2) If ”) P,(cos ,, Jy. .\(w,,) 
x[A, rei*¥+E, imi*] (50) 
and 


~ 2k(k—1)(k+ 1) —e 
2k+ 1 — Dy = 2we = >» 2, Pale Vir) w, P;,(cos ¢,,) 


K Sp 4 sail ¢,{A,, rei*+E,imi*]. (51) 
Eliminating D,_, between these equations we obtain another infinite set 
of equations 


x 


» 
> ¥a,(a%o*—V2 2)(57)? Pe. (cos 4, .;(t,)[ 4, rei* +E, im i*] — 0 
1 Bn Wh 
(k = 0,1,2,...) (52) 
if the field is to vanish at infinity. The remaining boundary condition is 


that the azimuthal resolute of the magnetic field shall vanish at the boun- 
dary. From (5) this implies that éw’/éz must vanish there and this gives 


7 2 Es —A, sin Baz +E,, cos Bn# A“ =) = 0 - 
2 on a a a 


n=1 Bn 
onr = a. By proceeding as before we obtain a third infinite set of equations 
for the constants, namely, 


wi “\ sind, Pi(cos ¢,,)J,.,(w, )[A, retv* +E, imi*] = 0 (54) 
w 


n 


(k = 1,2,...). 


Equations (43), (52), and (54) are the three sets of unknowns which are to 
be satisfied. 
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7. Numerical values 

In numerical illustrations of the orders of magnitude involved it will be 
useful to adopt values of the Alfvén velocity J, corresponding to the interior 
of the sun and the earth’s core. The values given by Alfvén (8) for the 
interior of the sun are of the order of 50-100 cm/sec, being based on a solar 
magnetic field of the order of 25 gauss. Babcock’s measurements show that 
this value is too large by a factor of order 10; hence we shall adopt here a 
value of } = 5em per sec. In the interior of the earth we may assume the 
earth’s field to correspond to that of a dipole field. According to Elsasser 
the radius of the core is of the order of 3-5 x 108 em so that the magnetic 
field at the equator of the core will be of the order of (6-4/3-5)* x 0-3 G or 
about 3G. Thus in the core itself the field will exceed this value and as a 
mean we shall adopt a value of 25 G. Taking also the density of the core to 
be about 8 gm_c.c. the corresponding Alfvén velocity Vi, is about 2-4 em sec. 
The observed equatorial velocity of the sun is 2 km sec, and assuming 
that the core rotates with the same angular speed as the earth as a whole 
the equatorial velocity of the core is about 0-25 km/sec. 


8. The effect of gravitational forces on the oscillations of a liquid 
globe 

We have already pointed out that the equilibrium of a liquid globe is 
neutral for internal displacements; hence gravity can only affect the period 
of the oscillations through the deformation of the surface, as is clear from 
equations (43) for the pressure boundary conditions. 

Now from (24) it follows that the order of magnitude of o will not greatly 
exceed Via since A, and y, are likely to be of order unity. Comparing the 
two terms in brackets in (43), remembering that 8,,, w,,, P,., P), are of order 
unity, the ratio of the second to the first is of order Gpa® V2. In the case of 
the sun we have p ~ 1, a = 6-9 10", VW, = 5 em/sec and it is abundantly 
clear that the gravitational terms in (43) are dominant. The same is true 
in the case of the earth’s core. The significance of this is that except when 
k = lwe must have e, = Oandifk = 1 we have arigid body displacement. 
Thus the surface of the globe remains spherical to a high degree of approxi- 
mation. The solution of the problem therefore also gives the oscillations of 
a liquid globe contained in a rigid fixed spherical boundary. 

The fact that the surface of the globe remains spherical enables us to 
effect a further simplification in (43) which may be written (omitting the 
factor k—1) in the form 


cn aii 2 
») ee VoBn Y=)! P*(c08$,,)Je.4(t,)[A,, te i* +E, im i*] = 0. 


(55) 


n Bn ; 2g 
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When & = | equation (43) reduces to (55) since 


P, (cos ¢,) B,, P\(cos $,, ) Ww, 


by (35) so that & takes all integral values in (55). 


9. The separation of the oscillations into even and odd modes 

As in the case of a non-rotating globe, so here the oscillations can be 
separated into even and odd modes. Equations (52), (54), and (55) are 
the three infinite sets of simultaneous equations in three infinite sets of 
unknowns. The condition for a self-consistent solution is the vanishing 
of an infinite determinant blocked in square matrices of order six. It is 
found that alternate blocks of 3 x 3 matrices in this determinant are zero, 
leading to two distinct sets of solutions in which w’ is either symmetric or 
antisymmetric with respect to the equatorial plane. The former solution 
will be said to correspond to an even mode; the other to an odd mode. The 
original infinite determinant now splits up into two vanishing infinite 
determinants, one for each mode, being blocked in groups of 3 x 3 matrices, 
Since the elements of the determinant depend implicitly on the period, 
these two vanishing determinants provide the corresponding period 
equations. 


The two determinants are derived by considering a pair of three infinite 


sets of linear equations for the constants A,, B,,, C,,, and for the constants 
E,,, F,, and G,,. 


” UJ 
Introducing the non-dimensional variable defined by 
a’ = aah, 
the cubic (23) becomes 


' 2)2/Q2 ) 2 " 

(o 2 B?)?(Bi, T Xn) fo . be 
where f = 2wy,a/N. (58) 
Using the values for w, and \, derived in section 7 for the sun and the earth’s 


core it is seen that f is a very large number in both cases. From (24) it 
follows that 


O<p,<o' <A,<@ and a, <»v, < o. (59) 


From (36)-(38) and (59) it then follows further that 
cos é., = —(o’?—2%)/fo’ 
C08 4, = (0p) /fo’ 
cosf, = (0’8+v5)/fo’ 
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and we have already noted that ¢,, is not real. Writing further 
A’, = «,(a?o'*— V2 A2)(27/z,)'A,, 
B,, x, (a?0'? — V2 42)(27/y,,)'B,, (61) 


Ch = a, (a?o"? + Vov7)(277/z,)'C, 


nr 


where z iz, and defining similarly E',, F;,, G, by replacing A,, B 


n? n? 
by E,, FP. G,, respectively in (61), we find that the equations for the even 


modes are 


y Py -2(COS é.. ) Joy . (2, \A . FY . P; 2k+ pa "n) die (Yn) Bi i 
ya) n@n Bn Jn 
nel 


Poy, 2(€08 F,, Lox (Zn )C, = 9, (62a) 
Pox, .2(€08 €,, + DY Pox .2(CO8 7) on 4(Yn) Bat 
n=l 
— FIP 5,908 f, )Tog (Zn) Cr 0, 


> Pr .1(COSE, oe (JAR LY Pax 1(CO8 9, ox .4(Yn) Bn t 
nel n=l 


2 | —)* Poy .1(€08 f,) ox .4(2,)C, = 0 
n 
whilst for the odd modes the equations are 
— (cos €, ) i — PP’ .(cos 7 a 
Dy eile) + > Bt aE ORY 
Aye, En Yn 


n=] n=l 


+>° - = Poy. .1(C08L,,)Tox .4(Zn) 


n=1 


(x, )E, T 2, Px + 1(COS8 nox ‘ (YF, 


¥ (=) Pag .1(008f,, log (2) @% = 0, (630) 
1 


n 


2x +2(C08 €,, oJ,,.,(2,,) £4 > Py. 9(€08 7, Joy. (4) F’,, 
n=l 


z ( )* Py .2(€08 Cn ons Zn) On 0, (63 ¢) 
1 


n 
where /,, is a modified Bessel function of order n. These equations hold for 
values of k = 0, 1, 2,... and the eliminant of the ratios of the constants from 
the sets (62 a—c) and (63 a—c) lead respectively to two infinite determinants 
to determine the free periods of the two modes. 
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10. Successive approximations to the fundamental periods 
(i) Odd modes 

The characteristic roots of the infinite determinant will be obtained by a 
method of successive approximation beginning with a 3 x 3 determinant. 
For the odd modes (63a, b, c) gives for the first approximation 


| A(x) Ayr Aly ey —Aledar | 
S21) Ayr) K(4) | =9. (64) 
J(x,)cosé, J(y,)cosn,  4,(z,) cosZ, 

This determinant factorizes but the only factors having zeros are J,(x,) and 

J;(y,). If 8 is the first zero ae then (60) gives 
‘2—NM)B = +fo'd, (65) 

and remembering that 2? = i +a?, we find that 
a 
and 8 > «,. In cases of geophysical and astrophysical interest we have 
already seen that fis large; thus one of these two roots is small and the other 

large, giving in fact 


= By(BP—ag) for fy (B?—ag)/B (67) 


whence, by (56) and (58), 


vi 2,/(p?— : 
= By(B?—at) 5 = or ee a) Wo. (68) 
Wo 


The first value of o agrees with the order of magnitude predicted by Cowling 
(1) for the fundamental mode. The other value of o is of order w,, that is, 
of the order of the angular velocity of the globe. But by (24) we have that o 
cannot much exceed V,/a, and if o is of order w,. this would imply that aw, 
must be less than Vj. In both the cases of the sun and the earth aw, is much 
larger than V. Thus the second value of o in (68) must be discarded. If, 
however, w,a<V,, that is the rotation is so small that hydromagnetic 
effects travel faster than the equatorial velocity, f is small and (66) gives 


o = r—oh(18), - 
thus both values of o’ are approximately equal in this case and the corre- 
sponding value of a is of the order of V,/a as found by Plumpton and Ferraro 
(2). But in this case we no longer get a continuous spectrum of frequencies 
as for the non-rotating globe. However, it is unlikely that such conditions 
will arise in nature as this would imply that the magnetic energy of the 
globe exceeds the energy of rotation. Taking f to be large we have that o’ 
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is given by the first of the two values in (67) and so is small. It then follows 
at once from (57) that, very nearly, 

2(62 +02) = f%o’?, (70) 
and that A,, 2%, = ¥, 2%, = fo’, gp, = Off -*). (7i) 
The last expression in (71) enables us to effect a further approximation for 
we have at once that (63a, b) reduce to the single equation 


n=1 n 


> — Pax (€08 Ean (nV Ent+ Y (—) Pap s1(08 Sn lox sa(2n)On = 9, (72) 
l 

and that B, = 0 very nearly. The equations governing the characteristic 

roots are then (72) and the equation 


’ cs 
SY Pox .2(cos €,, Igy 4(t, E+ > (—)* Poy. o(008 F,,) Lox. glZn)G, = 9, (73) 
n=1 
and in both equations k 0, 

The values of fo’ obtained by numerical calculation for successive 
approximations are given in the following Table 1. 


TABLE 1 




















It is seen at once that the convergence is not very rapid. This was to be 
expected since we are fitting a solution in cylindrical polar coordinates to a 
spherical boundary. From the table we should expect the fundamental 
value of fo’ to be eight or nine. 


(ii) Even modes 


Using the approximations (71), equations (62a—c) become 


> — Pox (eos &, )(7,)A,4 EI PS, (08 Lf, Lox (Zn )On = 9 (74) 
n=1 


and 


x 


> Pex-1(008 §,, Jo, (7, A, T > ( Poy .1(€08 flog (Zn )C) 0, (75) 
n n=l 


where, very nearly, Bi), = 0. The successive approximations to fo’ are 
given in Table 2. 


TABLE 2 





























OSCILLATIONS OF A ROTATING LIQUID SPHERE 299 


Thus the fundamental period for the even modes is about half that for the 
odd modes and we should expect the fundamental value of fo’ to be about 


21 or 22. 


11. Discussion 

The convergence of the above solutions is not very rapid, but since the 
6 « 6 determinant is in effect only the third approximation it has not proved 
possible to carry the computation further to the 8x8 approximation. 
The stream lines and lines of force have been drawn for each mode for this 
approximation. The figures will not be reproduced here, but it turns out 
that the coefficients of the exponential terms are a hundred times smaller 
than those of the trigonometric terms. This means that the exponential 
terms only take effect towards the poles of the sphere. Near the equator the 
functions are well behaved, indicating that the convergence of the series may 
be better near the equator. Since this problem was solved a digital com- 
puter has become available and the possibility of further computation is 
under consideration. It is hoped to report the results in a later paper. 

The periods of the fundamental even and odd modes for the earth’s core. 
taking fo’ — 9 and 21 respectively, V, = 2-4 cm/see and wya = 2-5 « 104 
em, sec, are 7 104 years and 3 x 104 years respectively. Thus the period of 
oscillation is comparable with the time of decay of the earth’s magnetic 
field, and so of little interest in connexion with the secular variation of the 
earth's field : this is most likely to be the result of dynamo action. 
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SUMMARY 
\ numerical method is deseribed for determining points on the neutral stability 
curve of a laminar jet. This involves the numerical solution of the Orr-Sommerfeld 
‘ (quation, yiven by 
d*d, dw i (d*d dd 
{u )| aa ) | 2a* v'h) 


4 «% 
dy? dy?” vh\dy'* dy?* 


with u ech®y and the boundary conditions 


lin 
verte 


The method used is to expand ¢ in a series of Chebyshev polynomials whose argu 
ment is the transformed independent variable ¢ tanh y, and solve the recurrence 
relations which exist between the coefficients in the expansion, 

An iterative technique is used for determining the parameters 2 and ree so that, 
for given a and zero imc, a non-trivial solution ¢ exists. A few points on the resulting 
eurve of R against a, the ‘neutral stability curve’, are found in order to estimate 
the ‘eritical Reynolds number’, R,, which is the minimum value of # for which such 


a solution exists 


1. Introduction 
THe Orr-Sommerfeld equation 


d*d Y dw i {[d‘d th 
( 4 c) "o® 2 2 ‘ 4 l 
u Fe b) ip? alae wate 4) (1) 


arises when the hydrodynamic stability of two-dimensional laminar shear 


flows is studied by investigating the growth of infinitesimally small 
disturbances. It is a fourth-order linear differential equation with variable 
complex coefficients for a complex function ¢ of a real independent variable 
y, w is a given real function of y, and «, R, and ¢ are parameters, « and FR 
being real and ¢ in general complex. In order to explain the significance 
of these quantities and to set the eigenvalue problem which is to be solved 
in connexion with (1), we now give a brief account of the approach to the 
+ Now at University of Adelaide, South Australia. 


Quart. Journ. Mech and Applied Math., Vol. XIII, Pt. 3, 1960) 
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hydrodynamic stability of nearly parallel flows leading to the Orr-Sommer- 
feld equation; for a more complete account see Lin (1, 2). 
Let x, y be cartesian coordinates, u, v the corresponding velocity com- 
ponents, and # the stream function, so that 
u = p,, v p,, (2) 
where suffixes denote partial derivatives. From the Navier-Stokes 
equations we have 
Vid + by, — PV, = vf, (3) 
where v is the kinematic viscosity. 
To investigate the stability for small disturbances of a basic steady flow 
given by ys — P%(x,y), we put 
fb = W(x, y) 4 p*(a, y, b), (4) 


substitute in (3) and neglect second-order terms in * and its derivatives. 


The terms corresponding to the main flow cancel out since they satisfy 


the steady Navier-Stokes equations and we are left with a linear homo 
geneous equation for p*: 
Vpt } ph) Vt yy” vy" | (VAye yp (V2? yp vV iy*. (5) 
We are concerned here with shear flows, for which the main flow is 
essentially parallel to one direction, say the a-direction; hence we now 
make the ‘parallel flow’ approximations. Great simplification accrues 
from this but the validity of the procedure has been questioned; however, 
the ultimate justification lies in experimental verification of the results 
obtained, and certainly for the case of boundary-layer flow along a flat 
plate agreement with experiment has been satisfactory. 
If we assume that the main flow is essentially parallel to the x-axis, we 
may firstly use the boundary-layer approximation and neglect the x 
derivative of any quantity connected with the main flow compared with 


its y derivative. Then (5) reduces to 


Vet hy ee Wy Pe = VAP". (6) 
Secondly, to the same order of approximation, we may ignore the x- 
dependence of the remaining main flow quantities and take their local 
values at some given value x, of 2. For most shear flows considered one can 
then use non-dimensional variables with respect to /ocal characteristic 
quantities, a velocity U(x,) and a length (2), so as to confine the depen 
dence on #, to a local Reynolds number R = U(x,)l(xq)/v. The equation 
for the non-dimensional * then becomes 


V2bt + w(y)V bt — w" (y)b® Aha (7) 
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where the non-dimensional ‘velocity profile’ w(y) is a function of y only, 
y now being the local non-dimensional coordinate perpendicular to the 
stream. The problem is thus reduced to the investigation of the stability 
of the ‘parallel’ flow with profile w(y) and Reynolds number R. Equation 
(7) is in fact identical with the equation relating to strictly parallel flow 
apart from the fact that for strictly parallel flow w"(y) is necessarily 
constant. 

Since the coefficients in (7) depend on y only, the other independent 
variables « and ¢ appearing only through derivatives, we may expect 
solutions of the form yt bly oie eb), (8) 
and assuming that any initial disturbance can be analysed as a Fourier 
series in 2 we may restrict attention to disturbances of this form. We 
consider only initial disturbances which remain small as 2 -» +, so that 
we further restrict a to be real. 

The Orr Sommerfeld equation (1) is obtained by substituting (8) in (7), 
and thus represents the condition that a disturbance of the form (8) shall 
be hvdrodynamically possible. Since the disturbance must vanish at 
infinity if the field of flow extends to infinity (or on the solid boundaries 
otherwise), 6 must also satisfy homogeneous boundary conditions of the 
form b $b 0 (9) 
at ¥ + & (or at finite values of y corresponding to the solid boundaries). 

[t will not be possible to find a non-trivial function ¢ satisfying all these 
conditions unless the parameters «, R, and ¢ satisfy a certain complex 


eigenvalue relation’, say F(x, R.c) 
- x, ,C€ 


For most physically reasonable protiles w(y) this relation appears to define 


e uniquely as a complex function of the real parameters « and R: 
e(a, R)+ie(a, R), 


and its imaginary part determines whether the disturbance (8) will grow 
(c; > 0) or decay (¢, < 0). Ife, — 0 for all values of «, the flow is stable 
with respect to infinitesimal disturbances; if ¢; > 0 for some value of x, 
we may say that the flow is unstable according to linear theory (see Lin 
(2)). 

For a given flow w(y) the (x, R) plane may be divided into regions where 
ce, < 0 and regions where c; > 0; these regions are separated by the curve 
given by ¢,(a, R) = 0, which is called the neutral stability curve. Lin (1) 
has shown that for several classes of flow this curve is in the form of a loop 
with its two branches asymptotic to lines of constant x as R-- & (see 
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Fig. 1), and that the inside of this loop corresponds to instability. Thus, 
fora given R > R,, where R, is the minimum value attained by R along 
the neutral stability curve, there is a band of wavelengths which are un- 
stable, but for R < R, all wavelengths are stable. Thus R, is the maxi- 
mum value of R for which all disturbances die out; it is the value of R 
at which instability sets in and is called the critical Reynolds number. Its 
determination is the main object of stability calculations. 


cj constant (<0) 
cj=0 


€;= constant (>0) 








In this paper we are concerned with solving this problem for the plane 
laminar jet with non-dimensional profile 
uw = sech*y. (10) 
For details of the non-dimensional form of Bickley’s boundary-layer solu- 
tion (3) see Curle (4), and for a discussion of the validity of the parallel 
flow approximations in this case see Tatsumi and Kakutani (5). Since the 
flow field extends to infinity in both directions, the boundary conditions 
to be satistied by ¢ are 


+ = 


lim 4 lim 4’ = 0. (11) 


Yori 
With the profile (10) and these boundary conditions the Orr-Sommerfeld 
boundary-value problem is completely symmetric in y, and consequently 
the odd and even parts of any solution ¢(y) are themselves solutions; they 
correspond to symmetrical (u even, v odd) and antisymmetrical (u odd, 
v even) disturbances, respectively. All previous work (see section 2) indi- 
cates that symmetrical disturbances give rise to a neutral stability curve 
which lies complete!y within the corresponding curve for antisymmetrical 
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disturbances and hence to a larger R,; this suggests that antisymmetrical 
disturbances would be the more common in practice, and this is borne 
out by observation. Indeed, this is to be expected, since the simplest 
mode of ¢ satisfying the boundary conditions is even, with a single maxi- 
mum at y = 0. Consequently only the case of even ¢ will be treated in 
detail. 

The numerical technique described in this paper (section 3) determines 
directly points on the neutral stability curve. We regard « and c; as fixed 
parameters (with c; = 0) and R and c¢, as two real eigenvalues to be 
determined, and by choosing a sequence of values of « we obtain a sequence 
of points («, R) on the neutral stability curve; R, may then be calculated 
by interpolation. Similar techniques may be devised for types of flow 
other than that of the laminar jet. The only other direct numerical treat- 
ment of the Orr-Sommerfeld equation has been that by Thomas (6) for 
the case w = 1—y*, which represents purely parallel flow under a constant 
pressure gradient between the parallel planes y = +1 (plane Poiseuille 
flow). Thomas replaced equation (1) and the associated boundary condi- 
tions by a set of finite-difference equations and, for given R and «, reduced 
them by Gaussian elimination for several values of the complex parameter 
c; the value of ¢ for which a non-trivial solution existed was then found 
by interpolation. This procedure was repeated with several different 
pivotal intervals and from the results the exact value of c was estimated. 
By this means a plot of ¢; in the (a, R) plane was obtained and then the 
critical Reynolds number determined by interpolation. 


2. Earlier work on the stability of a laminar jet 

The stability of the laminar jet has been investigated by several authors; 
the results of Curle (4, 7), Tatsumi and Kakutani (5), and Howard (8), 
are of particular interest. 


Using a method of approximation similar to one used earlier by McKoen 
(9), Curle (4) has calculated the neutral stability curve and the corre- 
sponding value of R, in the case when ¢ is even. He uses the known 
solutions 


x O,¢ 0.4 we sech?y: a 2,c= 


w = sech®y, 


of the inviscid equation (equation (1) with R=) and his method 
estimates the remainder of the curve for 0 <a < 2. The value of R, 


is found to be 5-5, corresponding to x = 0-22 and c = 0-05. 


~—- 
Tatsumi and Kakutani (5) obtain solutions of the Orr—Sommerfeld 
equation as power series in af for small values of aR, and inverse power 
series in aR for large aR. 
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The convergence of the series is too slow to permit the entire neutral 
stability curve to be evaluated by this means, but the branch of the curve 
for small « suffices to estimate the value of R, as 4-0, corresponding to 
a = 0-2 andc = 0-05. 

Howard (8) confirms that « tends to 0 or 2 as R -> 0, but shows that 
¢ = sech*y is not the corresponding solution at «a = 0. By employing a 
new perturbation procedure near « = 0, he estimates R, to be about 4. 

The case of odd ¢, corresponding to the problem of symmetric oscilla- 
tions, presents greater difficulty. Curle (7) observes that the inviscid 
equation has the solution, given by Savic and Murphy (10), 


x= l, c= a d = sech y tanh y 


but that no other simple odd solution has been found. He suggests that 
x = 0 again corresponds to a solution for infinite R, but that aR remains 
finite. Consequently, the right-hand side of (1) cannot be neglected when 
x» = 0, but may be approximated by —i(aR)~-'¢'¥. The equation thus 
modified was solved by one of us using a numerical method similar to that 
described in section 3 of this paper. It was found that the solution for 
realc gaveaR = 19, approximately. From other arguments, Curle showed 
that (1—a«?)R > 19, approximately, as a> 1. He then conjectured that 
the neutral stability curve might be approximated by a R(1—a?) = 19, 
which yields a minimum R of about 49 at a = 3-4. 

Subsequent computations on the lines of those described later in this 
paper for even ¢ indicate that the correct value of this minimum is some- 
what larger. However, since the results confirm that R, is considerably 
larger for odd ¢ than for even ¢, we do not pursue the former case further. 


3. Basis of the method 

In the present paper, a method for the numerical solution of ordinary 
differential equations (11) is used to solve the complete Orr-Sommerfeld 
equation in the jet case, with ¢ even. In outline, the method is as follows. 

After a suitable change of the independent variable, an expansion for 
¢é in Chebyshev polynomials is assumed, in which the coefficients are 
complex. Substituting this expansion into (1) yields a five-term recurrence 
relation between the coefficients in the expansion, with a, R, and ¢ as 
parameters. By an iterative scheme described in section 4, this relation is 
used to find approximate solutions of the equation, and hence to find the 
curve of neutral stability. The accuracy of the results depends on the 
number of terms included in the Chebyshev expansion of ¢. For a small 


number of terms, less than 15 say, the method is suitable for desk-machine 
5092 .51 x 





306 C. W. CLENSHAW AND D. ELLIOTT 


calculation, and in fact a first approximation to the neutral stability curve 
in the case of ¢ odd was found by this means, using 11 terms. In order 
to use many more terms, and to obtain enough accuracy to be able to 
extrapolate results corresponding to an infinite number of terms, the 
method was programmed for DEUCE, the digital computer at the National 
Physical Laboratory. 


4. Details of the method 

The independent variable is transformed by writing tf = tanhy. This 
has the effect firstly of changing the range of independent variable from 
(—x, +a) to (—1, +1), and secondly of transforming all the coefficients 
in the differential equation (1) into polynomials in the independent 
variable. This is a necessary preliminary to the application of the method 
(11) of expansion in the Chebyshev polynomials 7),(¢) = cos(n cost). 
A description of the properties and use of these polynomials has been 
given by Lanczos (12). 

For convenience we write C.{ f} for the coefficient of 7,(t) in the Cheby- 
shev expansion of a function f(t) for r > 0, and for twice this coefficient 
for r = 0. Let us suppose that 


so that d 4a, s a, T(t). (12) 
1 


n 


Let ¢ denote the sth derivative of 4 with respect to t. If we express the 
coefficients in the expansion of ¢ as 
Cid} = al, 


it can then be shown, as in (11), that for all r 


(13) 


and also that (4 ; (14) 


Using these results and the differential equation (1), we can find a five- 
term recurrence relation for the a, in the following manner. We have 

_ (dd) ; : 

( CA(1—#) 4b}. 

"\dy| rt } 


Equation (14) now gives 
CLEP} = f(a, + 2a +a ,), 


c 144 
"\dy| 


so that — }(a}) .— 2a +a{ 5). 
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Also equation (13) gives 
a) , = a+ 2(r—1)a,_, 
and at), = al) —2(r+1)a,,,. 


Substitution of these expressions in (15) yields 


fd 
Clr = (r+ 1)a,.,—(r—1)a,_,}. 


Similarly we obtain 


(9d) 


"\dy2) H(r+1)(r+-2)a, ,.—2r%a,+(r—1)(r—2)a,_2} 


and 


0 nA = pef(r+-1)(r+-2)(r+-3)(r+-4)a,.g—4(r+-1)(r + 2)(r?+-2r+-2)a,..+ 
. + 2r?(3r?+-5)a,—4(r— 1)(r—2)(r?—2r+2)a,_.4+ 
+(r—1)(r—2)(r—3)(r—4)a,_,}. 
Substituting these expressions in equation (1), and using the equation 
w = sech*’y = 1—?, 
we derive the required recurrence relation 
A(rja, 4+ Bir)a,_.+C(r)a,+ B(—r)a,,,+ A(—r)a,,4 = 9, (16) 
where 
A(r) = —K(r—1)(r—6)+ (r— 1)\(r—2)(r—3)(r—4), 
B(r) = 2K(r?—4r)+-4Ko?—B(r—1)(r—2)— 
—4i(r—1)(r—2)(r?—2r4+2+4-2 
C(r) = 2K(2—r?) + 2B(r?+ 2a2) + 2i(3r4+ Sr? 822+ Bat) 
and K = aR, B = 2K(2ce—1). 
The boundary conditions on ¢ are given by 


lim ¢@ = lim (1—#2)d = @; 
t-»+1 t--+1 


The first of these gives at once 


$4)—a,+a,—d,+... = 0) (19) 


hay+a,+a,+a,+... = 0) 


while the second is automatically satisfied if ¢ is finite; this is true if the 
Chebyshev series is convergent. The only remaining boundary condition 
to be applied is that the function ¢ is either an even or an odd function 
of y, and therefore also of t. In the even case a,,,, = 0 for all r, and in the 
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odd case aj, = 0. The conditions (19) can thus be reduced to the simpler 
forms 


$y +dg+a,+... = 0 (even) ) 


: (20) 
a,+a3+a,+...=0 (dodd) |} 

Equation (16) with r = 0, 1, 2,..., together with the appropriate member 
of equations (20), comprises an infinite set of linear algebraic equations 
for the a,. The problem is to find, for a given value of a, the values of K 
and § for which a solution is possible. 

It is known that ¢ is a continuous function of ¢ with bounded variation; 
hence its Chebyshev expansion converges. We may therefore assume that 
for some integer NV, the coefficients a,..,, 4, .9,..., are negligible to the work- 
ing accuracy. This reduces the equations to a finite set. 

We now consider a particular example in detail. The solution ¢ is taken 
to be an even function of t, and a, which lies in the interval (0, 2), is to be 
unity. We seek the real values of A, and B, which allow a non-trivial 
solution of equations (16), with r o.i.3 N, and the first of (20). 
First we solve equation (16) with r = N to find a,_,, taking arbitrary 
values for a,_, and a,, and neglecting a, for r > N. Equation (16) is 
then applied successively at r = N—2, N—4.,... to give ay _¢, ay_,,... until 
a, is obtained from (16) at r 4. We then have the following residual 
errors from the unsatisfied equations at r = 2 and r = 0, 


| A(2)+C(2)}a,+ B(2)a,+ B(—2)a,+ A(—2) 
and 2A(0)a,+2B(0)a,+C(O)ay = Ag, say. 


The boundary condition given by the first of (20) has also to be satisfied ; 

however, by summing the equations (16) from r = 0 to r = ©, it is found 

that 4a2(B+2K +4ia®)(Jao+a¢+0,+...) = pg+Ac. (21) 
20 2 4 2°*0 2 

Hence this condition is satisfied automatically if A, and A, are made zero, 

provided a sufficiently large value of NV has been chosen. 

A second trial solution is now obtained in a similar way, but starting 
from different arbitrary values of a,_, and ay. Any other solution of the 
recurrence relation formed in this way must be a linear combination of the 
two trial solutions, since there are only two independent solutions of (1) 
which remain finite as y > 0, that is, ast > 1. Denoting the trial solutions 
by superscripts (i) and (ii) we form the linear combination 


a, = ali 4 Bai), 
where the constant @ is determined by 


0 = AS? + 0rd). 
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-The remaining task is to find the real values of Ay and fy for which the 
(complex) value of 4 
I ) Ay = AL-+- BALD 


vanishes. 
We evaluate A, at three points in the (A, 8) plane and carry out a linear 
interpolation in the complex A,-plane. If we write, as an approximation, 


Ay = A+ BKy+CBy, 22 


then knowledge of A, at three (A, 8) points enables us to find the complex 
numbers 4, B, and C. The real and imaginary parts of the equation 

A+ BKy+CBy = 0 (23) 
then yield estimates of the required (A,8) which correspond to a zero 
value of Ay. These estimates are used to repeat the whole process until this 
requirement is satisfied. 

The corresponding set of coefficients a, represents the desired solution 
for the chosen value of V. If such solutions are obtained for various values 
of NV, an estimate may be made of the solution for NV = «, representing 
neutral stability. 

Finally, from the ultimate A and 8 we can compute 

R=oa IK and c= 4+ 4K. 24) 


5. Programme arrangement 

The programme for DEUCE corresponds closely to the scheme outlined 
above; only the following further remarks need be made. 

The parameters K and 8 may vary considerably in magnitude on and 
near the curve of neutral stability and we use instead the quantities ¢ and 
y = K~' which are more suitable for fixed-point arithmetic in a high-speed 
machine. The residual A, is now regarded as a linear function of ¢ and y 
in the neighbourhood of the trial solutions. 

The data required by the programme comprises the number N repre- 
senting the number of terms in the Chebyshev series, the quantity j = 320°, 
and initial estimates cy and y, say, of c and y. The programme then 
evaluates A, at three points, namely (¢9, yp), (Cg +4), Yo), and (¢g+-8,, ¥o+5e), 
where 6, and 6, are suitably chosen constants. From these three values 
of Ay a new (c,y) pair is calculated, which in general is better than the 
original estimates; that is, it will correspond to a smaller value of Aj. 
The output, from the programme, consists of the new (c,y) pair, the 
residual A, from the original estimate, and the corresponding value of ¢ 
at ¢= 1. As shown in section 4, this is zero only when N = o and the 
(c,y) pair is correct. The programme is so arranged that the output (c, y) 
can be used as a new input, and the cycle repeated so that yet another 
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improved result is obtained. The time taken to obtain one such pair from 
the previous estimate is approximately N/6 minutes on DEUCE, which 
has a multiplication time of 2 milliseconds. 


6. Numerical results 
Successive iterates for a typical point are shown in Table 1. Although 


TABLE 





18 





YN 





o'161 
o16g014 
"167879 


*295950 167982 


°) 
° 
296951 o'1679381 
° 
° 


*296950 


"167982 








these results are given here to 6 decimal places only, the arithmetic within 
the programme was carried out to double-length precision, corresponding 
to about 16 decimal places. Such a high working accuracy was necessary 
to offset the cancellation which may occur when two solutions are 
combined, 

Points on the curve for N = 24 were found at j = 1, 2, 18, 32, 72, and 
these form good first approximations for the corresponding points on the 
curve for VN = 48. Similarly, we obtain a curve for NV = 96, and in this 
case we take additional values for 7 = 4 and j = 8, since it is clear by 
now that the required minimum of FR occurs at a small value of j. The 
limiting values obtained for all these points are given to four figures in 
Table 2. 


TABLE 2 











cn YN Nn YN 





o1018 + 2°000 0°0322 + 1°499 
0°00 395 + 3°377 + 0°0554 17188 
or 1188 07392 
0*2009 o3gI5s 
o' 1680 + 0° 3004 +0°1634 0° 3017 0°1626 
0°07895 0°3935 0°07949 0°3933 0°07948 
0°01 840 + 0°5429 +o°o1 84! 0°5429 + O°O1 S41 




















Corresponding to each point in this table, we compute the values of Ry 
from the relations Ry = (ay)-! and « = (j/32)t. These are shown in Table 
3 against argument «. 
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TABLE 3 





a Ru 


= 
© 
> 





0°1768 
0°2500 
0°3536 
0" 5000 
0°7500 7°94 8:16 
1*0000 12°67 12°58 
1*5000 36°23 36°21 


NY OH WwW 
re) 














vu = @Gw 
2 -was 


“ 
_ 
- Oo 





We observe that for the larger values of « in the range considered, the 
values of cy and yy, and hence of Ry, vary only slightly with V. This 
implies that as a approaches 2, an accurate solution to the problem can 


a 








+ 
7"? 


Fic. 2 


be obtained with a quite small value of NV, as would indeed be expected 
if the solution of the equation (1) were tending towards the function 


w = sech*’y = 1—??. 


In contrast, the values of cy, yy, and Ry vary increasingly with N as a 
approaches zero, indicating that as « — 0 the solution of (1) is not a simple 
polynomial. This accords with the results of Howard (8) mentioned in 
section 2. 

Fig. 2 shows the curves of Ry against « for the three values of N; the 
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points labelled H are plotted from Howard’s approximation for large R, 
x = 0-954R-*+-5-4R-4+.... 


Tatsumi and Kakutani (5) also give the first term in this series, but 
with a different coefficient. Dr. J. T. Stuart, in a private communication, 
has pointed out their error, and agreement with Howard is now established. 


oz, 











We see from our curves that the required minimum of FR will exceed 
3°37; we may estimate it to be about 3-7. The corresponding value of « 
is about 0-25. 

In Fig. 3 are shown the corresponding curves of cy against a. Compari- 
son may be made with Fig. 2 of Tatsumi and Kakutani (5) which shows 
similar ill-behaviour for small «. 


7. Conclusion 

The present method gives results which agree reasonably well with those 
suggested earlier by Curle (4), Tatsumi and Kakutani (5), and Howard (8). 
We obtain a critical Reynolds number which lies between 3-37 and 4-0, 
probably about 3-7. The corresponding value of ¢ lies between 0-05 and 
0-10, and a is approximately 0-25. 


These results have been obtained by numerical solution of the complete 
Orr-Sommerfeld equation at a number of points on the curve of neutral 
stability in the neighbourhood of the required minimum. It is possible 
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that the equation may be solved for other types of flow by using similar 
expansions in Chebyshev polynomials, or other orthogonal functions. 
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SUMMARY 

The well-known result for the stability of a plane vortex sheet in an incom- 
pressible fluid is generalized to include magnetic fields and three-dimensional effects. 
The resulting current-vortex sheet is shown to be stabilized by sufficiently strong 
magnetic fields, and weaker fields are shown to give a partial stability. The effect 
of walls is examined, and in general, they are found to add to the stability of the 
flow. The dispersion relationships are found for several configurations of walls and 
current-vortex sheets, and conditions for stability are given. 


1. Introduction 
Tue stability of a sheet across which the tangential component of the 
magnetic field is discontinuous, has received a considerable amount of 
attention from astrophysicists and others. Usually this current sheet is 
assumed to have the shape of an infinite circular cylinder in its equilibrium 
condition. In this paper, the more elementary problem of a plane dis- 
continuity is considered, with the generalization that it is both a current 
and a vortex sheet. At such a sheet, discontinuities are allowed in the 
density and magnetic permeability of the fluid as well as in the tangential 
components of the magnetic field and velocity. Elsewhere, these quantities 
are assumed to be constant, and the streamlines and lines of force are 
straight lines lying parallel to the surfaces of discontinuity. Streamlines 
and lines of force which are separated by a surface of discontinuity need not 
be parallel, corresponding to the presence of both axial and azimuthal fields 
in the cylindrical case. It is assumed that the fluid is incompressible, and 
that the effects of viscosity and of finite electrical conductivity can be 
neglected. 

The equations to be satisfied when there are small departures from this 
equilibrium situation are derived, and the general solutions found for the 
case in which the current-vortex sheets are given a small sinusoidal disturb- 


ance. The conditions with which these solutions must comply at a current- 
vortex sheet, and at a wall, are given in a suitable form, and used to find the 
dispersion relationships in several examples with various combinations of 


+ The author is a member of the New Zealand Scientific Defence Corps. 
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current-vortex sheets and walls. The simple geometry assumed for the 
equilibrium state of the field and flow, allows quite complicated solutions 
to be obtained with little difficulty, and there is the advantage that a com- 
parison can be made with results already well known in ordinary hydro- 
dynamics (Lamb (1), section 232), where the surfaces of discontinuity are 
vortex sheets only. These results are given by Lamb for two-dimensional 
cases. The generalization to three dimensions, together with the effect of 
walls, can be obtained from the results of this paper, if the magnetic terms 
are ignored. Some particular cases of the results given here have been 
obtained by Michael (2). 


2. Linearization of equations 
The magnetohydrodynamic equations for an incompressible fluid, with 
terms involving the viscosity and magnetic diffusivity omitted, aret 


6 4(V.0)B ~ (B.Y)V, 
; 


V.B= 0, 
ov l B*\ 1 
—+(V.V)V = — v(P+ = )+ (B.V)B, 
et p “H} ep 
V.V=0, 
where B is the magnetic induction, V the fluid velocity, P the pressure, 
and yu and p the permeability and density of the fluid respectively. Taking 
a system of cartesian coordinates (x, y, z) such that current-vortex sheets 
and walls are planes of constant z, the field, velocity, and pressure can be 


written as B — B+b = (4, B,0)+(a,b,¢) 
V = Vav = (U,V,0)+(u, v, w) }- (5) 
P= P+p 

The equilibrium conditions are given by the first term in these expressions, 
and the second term represents a perturbation caused by a disturbance to 
the current-vortex sheets. Since A, B, U, V, and P are assumed to be con- 

stant, the linearized forms of the above equations are therefore 
(5+ Us + v <p = (4 c+BSw, (6) 
V.b = 0, (7) 
at ta =P tet a 
V.v=0. (9) 


+ M.K.S. units are used throughout. 
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It will be assumed that the lines of force and the streamlines coincide so 
that there is no electric field, and the constants A, B, U,V are not inde- 
pendent. In fact BOY 


4 Uv A, -_ 


and tan~'A is the angle made by the streamlines and lines of force with the 
x-axis, This condition can always be brought about by a suitable choice of 
moving axes; for the electric field, E, satisfies 


E+VAB = 0, 


and if axes are chosen to make E zero in one region, then the condition that 
the tangential component of E should be continuous at a surface of dis- 
continuity ensures that E is zero everywhere. 


3. Solutions for sinusoidal disturbances 

Consider the case in which the z-coordinate of a current-vortex sheet is 
disturbed by an amount 

8O 806 ¢ Kottke+ly) (11) 

where k, / are arbitrary. The amplitude 5, must be small compared with 
any typical length in the equilibrium configuration. The quantity s, which 
is of unit magnitude, is inserted to provide differences of phase where 
required, Since the equations derived in the last section are linear, solutions 
can be superimposed, and so it is sufficient to consider a sinusoidal disturb- 
ance, as any arbitrary disturbance can be analysed by a Fourier series 
into a sum of components of this type. 


Solutions are sought having the form 
b = Kéde'™ 
v = Xde™ (12) 
p pDbe'™ 
where K and X are constant vectors such that K = (AK, L, M), 


X = (X,Y, Z), and ) isa constant. It is easily found by substitution into 
equations (6), (7), (8), and (9), that 
k?+-[?4-m? = 0, (13) 
BD 
(wa)? 
' D ‘ 
X k, (15) 


(cw + a) 


k, (14) 


where 


kU+lV U(k-+-Al), (16) 
kA+1B = AalU, (17) 





THE STABILITY OF PLANE CURRENT-VORTEX SHEETS 317 
and k is the vector (k, 1, m). Equation (13) gives two values for m, namely, 
m == +4(k®+-[?)t = +in. (18) 
Hence, in any region of the fluid which does not contain a discontinuity, 
there are two linearly independent solutions for the disturbances, and the 
general solution is 
b = [(k, 1, in)D’e-"*+-(k, 1, —in) D"e"* |B3 ~ (w+ a)’, (19) 
v = [(k, 1, in)D’e-™ +-(k, l, —in) De" |b (w+-a), (20) 
p = [ D'e-"=4 D"e"=|p8, (21) 
where D’, D” are arbitrary constants. 
The solutions given by equations (13), (14), and (15) are not the only 
oscillatory ones; there are also the Alfvén wave solutions, such that 
pp(w+a)? = B?, 
and ppD+B = 0. 
These solutions are not relevant to the present work however, as they 


correspond only to real values of w, whereas complex values of w are required 
if any instability is to occur. 


4. Boundary conditions 
If the fluid has a solid boundary in the planez = h, the normal component 
of the velocity must vanish in this plane. Hence, from equation (20), the 
boundary condition is seen to be 
D'e-mh Der, (22) 
This implies that c is zero at the wall, and so no further condition is needed. 
The conditions which must be satisfied at a current-vortex sheet are 


(i) the total pressure, (P+-B?/2u), must be continuous, 
(ii) the normal velocity on each side must be equal to the normal velocity 
of the sheet, 


(iii) the normal component of the magnetic induction must be con- 
tinuous. 


The surface of discontinuity is taken to be 
z= h+s86. 


Hence the linearized forms of the above conditions are 


( _ Aa+ =) ( Aa-+ *) 
p- | = ([p+——— onz =/h, 
ad 1 B 2 


e (* 1,2, 7% 


i... 1 "ac 
OX cy 


on z =A, 





. AXFORD 
, 6 
ox 


08 08 0d 
] ( 8A 8B on z h. 
1 2 


oY Ox oy 


Ho) on z =A, (26) 
cy 


The first three conditions can be written, using the forms (19), (20), and 
21), as : 
2 2 
Py Bi . [Dye nh 4 Dien") Pe Pa - | Dye nh 4 Dye"), 
fy (w 4 x)? fig(w 4 xy)? 
28) 
y 


~(w+a,)? = Die-™— Di e™, (29) 
n 


fe vy)? Dye nh Dye. (30) 
n 


Using the third component of equation (1), and the conditions (25) and (26), 
it can be shown that both sides of (27) are zero. Thus the condition on the 
normal component of the magnetic field is satisfied automatically if 
(25) and (26) are satisfied. 

For cases in which the fluid extends to infinity, the condition which must 
be satisfied is that the disturbance must vanish at large distances from its 
source. Hence, since n is positive, 

DD” — Oif the fluid extends to +00 ) - 
| (31) 
p 0 if the fluid extends to ae) 


5. The plane current-vortex sheet in an infinite region of fluid 


Consider a current-vortex sheet in the plane z — 0, with fluid extending 
to infinity on either side. Let subscripts | and 2 refer to the fluid in z —- 0 
and z < 0 respectively. The conditions (31) give 


A Dy, 0, (32) 


and the remaining conditions, (28), (29), and (30), become 


2 


22 2 
Py Py Dy, |. B3 LY, (33) 


fylw- x)? pa(w ' x, )? 


8 


~ (w+ a4)? I, (34) 


n 


* (w+ ag)? Di; (35) 
n 
The phase term, s, has been retained here, although it has no significance 


in cases with only one current-vortex sheet. The constants 1), and Dj 
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can be eliminated from these equations, leading to an equation involving w, 
which is the required dispersion relation, 


(py + Pa)u® + 2(ay py + gpg) + |. at +p, 0g) — (A Fi)| = 0. (36) 
Pi Fe 


It is convenient at this stage to introduce the dimensionless quantity 2?, 
defined as Q2 — A*/ppU*? = BY/ppat. (37) 


This quantity, expressing the ratio of the magnetic and kinetic energies 
per unit volume, or alternatively, the ratio of the Alfvén velocity to the 
fluid velocity, is analogous to the Mach number in ordinary compressible 
flow. The dispersion relationship (36) can now be written 


(py + p_)w® + 2(ay py + ag pg)w +p, (1-22 )a? +-p,(1 —Q3) a9 0. (38) 


For stability, it is necessary that for all modes, the imaginary part of w 
should not be negative. Hence the discriminant of equation (38) must be 
positive, and therefore 

(Py + Pa) pr 82 of + py 823 0G) > py polory— 1g)? (39) 
This can be regarded as defining, for particular values of the parameters 
involved, directions of the disturbance vector k for which there is stability. 


The sufficient conditions that there should be stability for all directions of 
the disturbance are 

(pr +p2)Q2i > po, (py + p2)823 > pr, (40) 
and PrP2 < [ (Pr +p2)62i — poll (r+ p2)823 —py|- (41) 
The regions in the (Q?, Q3) plane, in which there is stability for all disturbance 
directions, are shown in Fig. | for various values of the density ratio p,/p,. 
Note that if Q? > 1 on both sides of the current vortex sheet, there is 
stability for all disturbances, but if Q? < 1 on both sides, there is some 
instability. This could be expected, for the magnetic fields produce all the 
stabilizing effects, and if Q? < 1 on one side of the current-vortex sheet 
(equivalent to ‘supersonic’ conditions), the field in this region cannot 
adjust itself rapidly enough to prevent the growth of small disturbances. 
Thus it is necessary that Q? > 1 (‘subsonic’) on one side of the current- 
vortex sheet at least, although this is not a sufficient condition, as can be 
seen from Fig. 1. 

The results for the ordinary hydrodynamical problem can be obtained 
directly from the above by putting Q, = Q, = 0. In particular, (39) shows 
that a vortex sheet is always unstable in the absence of a magnetic field. 
With weak magnetic fields, partial stabilization is produced, for there will 
be stability for a small range of directions of k. As the magnetic fields are 
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increased in strength, this range becomes larger, and eventually all direc- 
tions are stable when the conditions (40) and (41) are fulfilled. In order to 
show this behaviour, it is convenient to make a rotation and change of 
velocity of the frame of reference so that 
A, = —A, = tané, 
U, = —U, = U. 


Seqteeeerecs 
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Stable for all pi /P2 


a RD Re A a 
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Fic. 1, Regions of stability in the (Q7,Q5) plane for a single current-vortex 
sheet with arbitrary orientation of the field and flow on either side. Each 
curve represents the stability boundary for a certain value of the density 
ratio: (4) py/p, ~ 9, (b) py p, = 4+ (©) Pipa 1, (d) py/p2 = 3, (€) pip, = ®. 


The hatching denotes the unstable side of the boundary in each case. 
Thus if M—"TPspe = yn — fit Paps 
Pe Pi 
and the disturbance vector k makes an angle ¢ with the z-axis so that 
tand = I/k, equation (39) becomes 
(tan @ tan d)?(M+-N)—2(tan étand)(M—N)+(M+N—4)> 0. (42) 
If # is zero, the condition for stability is simply M+ N > 4, which can be 
written ‘ 
13 ‘ | _ PrP2 ye 


= Fo = (43) 
— 6s Pit Pe 


= 
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This is the result given by Michael (2) if p; = p, and p, = py. There is 
stability for all values of tang when @ = 42 provided (M-+-N) does not 
vanish, for in this case the vortex sheet disappears. Fig. 2 shows the 
condition (42) in terms of tand and M, for various values of tané. Two 
cases are shown, one with N zero, and the other with M and N equal. Note 
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Fic, 2. Two cases of the condition given in equation (42) for the stability 
of asingle current-vortex sheet, one with N = O and the other with M = N. 
In each case the regions below the curves represent unstable conditions, 


that the stability is better in the second case, for if N is zero and 
0) < 6 < 4m there is always one direction given by tang = cot @, in which 
there is instability, however large M may be. On the other hand, if 
M — N > 2 there is stability for all directions of the disturbance, In 
general there is instability if tan@tan¢ lies between two limits, which 
differ by 4M 4+N—MN)}/(M+N). 


6. A single current-vortex sheet in a bounded region of fluid 
In this case the current-vortex sheet lies in the plane z = 0, and there 
are walls at z = h,, andz = —h,. If subscripts | and 2 are used to denote 
the fluid in z > 0 and z < 0 respectively, the dispersion relationship is 
found to be 
py| (w+ a)? —QF af | coth nh, + pol (w+ a,)? —Q3 ag] coth nh, = 0. (44) 


This is the same as (36), but with p, replaced by p, coth nh, and p, replaced 
5092.51 Y 





322 W. I. AXFORD 


by p,cothnh,. The condition for stability is therefore found by making 
these substitutions in (39). Note that (36) and (44) are identical when 
h, > 2%, h, > ®, and also when h, = A,; hence it can be argued that sym- 
metrically placed walls have no effect on the stability. 

The dispersion relationship for the case in which only one region (say 
region 1) is bounded, can be found from (44) by letting , + «, and the 
stability condition is found from (39) by substituting p, coth nh, for p,. For 
small values of nh,, the condition for stability is approximately 

pt PY 
MiP: Pe2Pe 
and hence, provided 8, ~ 0, there is stability for disturbances of sufficiently 
small wave number, Thus the effect of a wall is to stabilize the current 


BF * Py pynh, (a4)? 


(45) 


vortex sheet against disturbances with wavelength O(/,) and greater. 
If there is no magnetic field in the bounded region of fluid, there is in 
stability, but it is not as severe as it would be with no wall. 

The configuration with two walls is not as stable as that with only one, 
and if n(h, + hy) 1, there is no relation similar to (45) indicating stability 
at sufficiently small wave numbers. However, if nh, » 1 and nh, l, 
there is stability in a certain range which does not include the smallest 
wave numbers. 


7. Configurations with two current-vortex sheets 


Consider the symmetrical equilibrium situation in which a region 2 is 
separated from regions | and 3 by current-vortex sheets at z thy. 
The field, velocity, and fluid properties are assumed to be the same in 
regions | and 3, which are bounded by walls at z thy (hy > hy). There 
are two types of overall disturbance, one being symmetrical (with the 
disturbances in the current-vortex sheets exactly out of phase), and the 
other asymmetrical (disturbances in phase). 

The solution for the symmetrical case is found by taking s loon 

h,, and s lon 2 h,. The resulting dispersion relationship is 
pil (wm u,)* (2? i tanh nh, + pol (w ; v,)? QQ x3| tanh n(h, h,) 0. (46) 
For the asymmetrical case, 8 is taken to be +1 on z +-h,, and the dis 
persion relationship is then found to be 


pil (w- x,)? Q)? yi| coth nh, + pol (w 4 x»)? (22 y3| tanh n(h, h,) 0, (47) 


These results agree with those given by Lamb when there is no magnetic 
field. It is a necessary condition for stability that the discriminants of 
equations (46) and (47) (regarded as quadratics in w) should be positive. 
These conditions are found by substituting p,tanhn(h,--h,) for p,, and 
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either p, tanh nh, or p, coth nh, as the case may be, for p, in (39). Stability 
diagrams such as shown in Fig. 2 are rather complicated in this case but 
portions of the curves for which (k*-+-/*)* is either very small or very large, 
can be found without much difficulty if suitable approximations are sub- 
stituted for the functions tanh and coth. 

For the case in which a, = B, = 0, representing a jet surrounded by a 
magnetic field, the configuration is stable (since tanh nh, < coth nh,) if 


2 
tanh nh, > Pal Ha 1 tanh n(h,—h,), (48) 
Pi 1 
2 


and 8, # 0. Thus a sufficient condition is that of < Pi/y,p,. If nh, < 
then (48) is approximately 


2 
prhy > pylhy—hj) a - | 
By 
but if h, + «© and nh, < 1, the approximate form is 


2 
xX, 
np, hy > py a |. (50) 
1 
Hence if 44, p,; «3 > Bj, all wave numbers less than a value given by (50) are 
unstable in the absence of walls, but with walls there is stability for wave 
numbers < 1/h, provided 


hy « iy} : pi : mt (51) 


Poh (Ha Py x3 /B7) 

In more general cases with two current-vortex sheets, where the equi- 
librium configuration is not symmetrical, the disturbances to the current- 
vortex sheets must be allowed to have different amplitudes. The ratio of 
these amplitudes, as well as the dispersion relationship is found from the 


conditions of section 4. 


8. Configurations with parallel streams 

It has been observed in sections 6 and 7 that the dispersion relationships 
(44), (46), and (47) are of the same form as (39) but with p, and p, replaced 
by terms involving the wave number n, say p,(m) and p,(n). A general 
result covering all these cases can be obtained when the streams are parallel 
in the undisturbed state (i.e. A, — A,); for the factor (4+ Al)? may then be 
removed from the general condition for stability which reduces to 


[p,(m) +-py(n) |l py(m)Q? VF + py(njQ3 U3] > py(n)pyin)(U,—0,)%. (52) 

This condition is obviously satisfied if U U,. When 1 4 U,, it is con- 
venient to choose a frame of reference such that ( -U, = U. This 
can be done without affecting the electric field since V ~ B is unaltered 
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by the translation. Equation (52) is now a quadratic form in p,(n) and p(n) 
which must be positive definite if there is to be stability for all n. The 
required condition is 

(02? — 023)? < 8(Q2 +03) — 16, (53) 


and the region of stability in the (QQ7, 3) plane is therefore the interior of a 


parabola with its vertex at Q? — Q3 — 1. The axis of the parabola is the 


line Q? — Q2, and the curve touches the coordinate axes at the points 
(0,4) and (4,0), 
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SUMMARY 
The theory is developed of oscillatory flow of the most general linear eclastico- 
viscous liquid in a concentric-sphere elastoviscometer, in which the outer of two 
concentric spheres undergoes forced harmonic angular oscillations about a vertical 
axis and the inner sphere is constrained by a torsion wire, [tis concluded that such 
an instrument would distinguish different clastico-viscous liquids in much the same 


way as the coaxial-cylinder elastoviscometer, 


1. Introduction 

In the present paper the motion is studied of an elastico-viscous liquid 
contained between concentric spheres. The outer sphere is assumed to 
undergo forced harmonic angular oscillations about a vertical axis through 
its centre, while the inner sphere is constrained by a torsion wire. With 
this apparatus it should be possible to measure the ratio of the angular 
amplitudes of the two spheres and also the phase-lag between the spheres. 
End effects would not need to be considered, and from this standpoint the 
concentric-sphere elastoviscometer would have a theoretical advantage 
over the corresponding coaxial-cylinder instrument. A further advantage 
arises from the fact that the general equation governing motion in the 
concentric-sphere instrument can be expressed in a form suitable for com- 
putation without introducing any approximation into the theory. This is in 
contrast to the corresponding theory for the coaxial-cylinder instruments 
(1, 2). An important practical advantage of the spherical over the cylin- 
drical instrument might well be that smaller samples could be used; and 
it is suggested that this might make it worth while to try to overcome the 
obvious constructional difficulties if the instrument is shown to be useful 
in principle. 

Although a number of authors have considered the use of oscillating 
spheres in rheological investigations (3-6), no concentric-sphere elasto- 
viscometer has yet been described in the literature. The mathematical 
analysis presented below, which is very similar to that given elsewhere 
for the coaxial-cylinder elastoviscometer (1, 2), is therefore undertaken 
to assess the potential value of a spherical instrument in rheological in- 
vestigations. It is shown that measurements of amplitude ratio alone 

[Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 3, 1960] 
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in a spherical instrument would not be sufficient to resolve the ambiguities 
of interpretation arising from comparable measurements in a coaxial- 
cylinder elastoviscometer (2), and it is concluded that the concentric- 
sphere elastoviscometer would distinguish different elastico-viscous liquids 
in very much the same way as the coaxial-cylinder elastoviscometer. 


2. Flow between concentric spheres 

The general linear equation of state for an elastico-viscous liquid has 
been discussed in the preceding paper (2). We can consider the stress at 
any point as the superposition of two stress systems, one (p;,,) associated 
with distortion and the other (— pg;,) with dilatation, and we write 


Pik = Pik—PVix (1) 
where p is a scalar and g,, is the metric tensor. The general linear elastico- 
viscous liquid has an equation of state of the form 


t 
Pix = 2 | F(t—t'eR(t’) dt’, 


£ 


x 
e 


| N(r) | 


where Y(t—t’) tt) dr, 


0 


being the relaxation function and N(r) the distribution function of 
relaxation times. If the liquid is incompressible, there is the added condi- 
tion that eft) 0 (all p). (4) 

We desire to investigate the motion when a general liquid of this type 
is contained between concentric spheres and subjected to a motion brought 
about by submitting the outer sphere to forced harmonic angular oscilla- 
tions about a vertical axis through its centre, the inner sphere being 
constrained by a torsion wire. 

We shall express all tensor quantities in terms of their physical com- 
ponents in spherical polar coordinates (r, x,¢), d being the azimuth angle 
about a polar axis drawn vertically upwards. In terms of these coordi- 
nates, two-dimensional axially symmetric motion is assumed in which the 
velocity components are [0,0,rsin y(éw/ér)| so that w(r,t) is the local 
angular velocity about the vertical diameter. The only important equation 
of state is that involving rd’, which becomes 


t 
~ Cw 


rp’ = rsin x | Y(t—t’) — dt’. 
ér 


x 
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This equation has to be solved in conjunction with the equation of motion 
ord’ 3, : Cw : 
= + 1b = preinx— (6) 
where p is the fluid density, and the appropriate boundary conditions. 
If n is the frequency of oscillation, a solution is possible of the form 
rp’ = re f(r)sin ye?7"™, w = re F(r)e27i", (7) 


where f(r) and F(r) incorporate appropriate phase factors. Substituting 
in (5) and (6) and eliminating F(r), we find 


2 
(n+ me 1. rt—6)f 0, (8) 


— 2minp —2rinp 
x - a * (9) 


[ W(gje*7iné dg | {N(r) de/(1+-2ninz)} 


0 0 


where 


The solution of equation (8) can be expressed in terms of Bessel func- 
tions of complex argument in the form 


f =r {aJ,(ar)+bJ_,(ar))], (10) 


where a and + are constants to be determined from the boundary condi- 

tions. The corresponding expression for F(r) is 

’ mx 
F = * _[aJ,(ar)—bJ_,(ar)}. (11) 
2rinpri 
It is convenient to express the boundary conditions in terms of the 
angular amplitudes of the two spheres. If 6, and @, are the amplitudes 
of angular oscillation of the inner and outer spheres (radii r, and r, respec- 
tively), then 

F(r,) = 27inO,, F(r,) = 27in8,e", (12) 


where c is a real constant defining the phase-lag between the spheres. A 
small positive c implies that the inner sphere lags slightly behind the 
outer. 

We also have as a boundary condition the equation of motion of the 
inner sphere, which reduces to 


Sori f(r,) = (K—4n°n?1)6,, (13) 


where / is the moment of inertia of the inner sphere about its vertical 
diameter and K is the restoring constant of the torsion wire. 
The solution of equation (8), subject to the boundary conditions (12) 
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and (13), is most conveniently expressed in the form 

al J_,(arg)J, (ar, )+J (ory Jj (ocr) | + ui] J; ( xr) y(27rg)— Jy(ary)J. (or, ) 
ae (r9/r3)'[ J_y(ar,)J,(ar,)4 Jj(ar,)J_,(ar,)] 

(14) 
3 K 
where S = A; —I}. (15) 
” 8apri 4n*n? 

This equation permits calculation of both the amplitude ratio # (= @,/6,) 
and the phase-lag c. Using the exact expressions for the Bessel functions 
of orders equal to half an odd integer, in terms of algebraic and trigono- 
metrical functions (7), we can recast equation (14) in a form more suitable 
for yy 


rm — b tes 
(") || {(: r;) j —_ j x-+ S sin(r,—?r,)x-+ 
11s “TL! 1 7.% 


3(r2—r,) S(re—r 
2,2 . we aaa 
TTT s "i" 


| cos(r, r,)x|- (16) 


Unlike the corresponding equation for the coaxial-cylinder instruments 

(equation (27) of (1)) which is valid only for large arguments of the Bessel 

functions, equation (16) is exact and can be used for all values of «. 
For very low frequencies 

ee (rB ry)Ki (17) 

de ° l6r*rird nen 

where np is the limiting viseosity at small rates of shear (cf. 2). This 


implies that 
implies th Pe her as n—>O, (18) 


a result which is independent of the relaxation spectrum and is true for 
all linear elastico-viscous liquids. 


3. Theoretical predictions 

Since no actual concentric-sphere elastoviscometer is available for com- 
parison, numerical calculations must be made by assigning hypothetical 
values to the physical constants of the apparatus. It seems reasonable 
to take values comparable to those used in the cylindrical apparatus of 
Oldroyd, Strawbridge, and Toms (8); for in this way it is also possible 
to compare the sensitivity of a concentric-sphere apparatus with that of 
the corresponding coaxial-cylinder instrument. Values of A and J are 
accordingly chosen to correspond to the A (cylinder II, torsion wire 8) 
set-up of Oldroyd, Strawbridge, and Toms (8); i.e. K = 2-79 x 10® dyne 
em rad-!, and J = 2468 gm cm?.__r, is chosen to be 4 cm and two radial 
gaps are considered with r,—r, taking the values 0-1 cm and 0-2 em. 

In order to make specific theoretical predictions, we consider first the 
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Fic. 1. The variation of 3(n) with B when p 
(r2—T; 0-O1 em.) 


prasemapiunieanens = 
n(e/s) 20 





Fic. 2. The variation of ¢e(n) with B when p 
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two-constant prototype liquid introduced in (2), with a relaxation spec- 
trum of the form 


N(r) 





TJo 
p 


For this simple prototype 











47*n*pB 


— (20) 
no log(1 + 27inB) 


Figs. | and 2 contain some (#,n) and (c,n) curves for this prototype. 
The curves show the same general features as those relating to the coaxial- 
cylinder instrument (Figs. 1 and 14 of (2)), each (#,n) curve showing a 
peak at a frequency greater than the natural frequency », with # passing 
through unity at a value of n almost exactly equal to ny and each (c, 7) 
curve passing through zero also at a frequency almost exactly equal to np. 

We conclude that the concentric-sphere elastoviscometer would distin- 
guish different elastico-viscous liquids in much the same way as the 
coaxial-cylinder elastoviscometer. This conclusion is reinforced when con- 
sideration is given to the possible use of this elastoviscometer in resolving 
the ambiguities of interpretation discussed in (2). There it was shown 


that experimental (3,) curves for certain dilute polymer solutions could 


be interpreted in more than one way—in terms of a discrete relaxation 
spectrum or in terms of a continuous relaxation spectrum. It was also 
shown that measurements of phase-lag at high frequencies might in 
principle partly resolve such ambiguities. It is, therefore, of some interest 
to plot (#,n) and (c,n) curves for the concentric-sphere elastoviscometer, 
corresponding to two spectra—one continuous and one discrete—which 
would be indistinguishable in the amplitude-ratio experiments of Oldroyd, 
Strawbridge, and Toms. In this way, it is possible to determine whether 
ambiguities can more readily be avoided by using a concentric-sphere 
instead of a coaxial-cylinder apparatus. 

We focus attention on the particular spectra that played a prominent 
part in the discussion of (2). The discrete spectrum is represented by the 
three constantst y, = 7-9, A, — 0-065, A, = 0-015; and the continuous 
spectrum by the three constants ny = 7-9, 0 = 0-13, B = 0-18. Figs. 3-6 
contain the theoretical (3,n) and (c,n) curves for these spectra in two 
hypothetical experimental arrangements. It is clear that amplitude-ratio 


The notation is that of (1). 
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Fic. 3. The relation between & and n in the case of both continuous (full 
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Fic. 4. The relation between ¢ and n in the case of both continuous (full 
line) and discrete (broken line) spectra. (r,—7, 0-1 cm.) 
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experiments would not be able to distinguish the two spectra, but that 
measurements of phase lag at high frequencies, taken in conjunction with 
them, would be able to distinguish them. An identical state of affairs 
exists in coaxial-cylinder experiments and the above conclusion is con- 
firmed. 
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SUMMARY 


The stresses and displacements within an elastic-plastic curved bar of perfectly 
plastic, non-hardening, compressible material subjected to a pure bending moment 
are calculated for both the Tresca and von Mises yield conditions and their associated 
flow rules. Plane-stress conditions are assumed and the applied bending moment 
is determined in the four cases which occur for either yield condition, Some numerical 
results are presented and the results for the two yield conditions are compared. 


1. Introduction 

SHAFFER and House (1, 2, 3) have given the solution to the problem of 
the elastic-plastic bending of a wide curved bar of incompressible material 
under conditions of plane strain. The solution to this problem when the 
material is compressible has been presented by Eason (4). The corre- 
sponding problem under conditions of plane stress was investigated by 
Shepherd and Gaydon (5). They considered a fully plastic curved bar of 
rigid-plastic material for both the Tresca and von Mises yield conditions. 

The problem considered here is that of the elastic-plastic bending of a 
curved bar by end couples under plane-stress conditions. Solutions are 
obtained for both the Tresea and von Mises yield conditions and their 
associated flow rules. Compressible perfectly. plastic material which does 
not work harden is assumed. Conditions are envisaged in which the bar 
is bent in its own plane and instability is assumed not to occur. The 
assumptions are made throughout that the plate thickness does not change 
appreciably. Whilst the plastie-strain components remain of the same 
order of magnitude as the elastic-strain components this assumption is 
unlikely to introduce serious errors. 

Four possible cases are found to occur for both yield conditions, these 
being distinguished by the number of elastic and plastic regions present 
in the plate. Expressions for the stresses, displacements, and applied 
bending moment are presented in all the possible cases, and some numerical 
results are given. Comparison is made between the two yield conditions 
for the same applied bending moment. It is found that for cases in which 


the stresses do not differ signiticantly there may be quite large differences 


in the area of the plate which is plastic. Two integrals which are of impor- 


(Quart. Journ, Mech. and Applied Math., Vol. XIII, Pt. 3, 1960) 
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tance in the solution for the von Mises yield condition are tabulated. The 
numerical work was programmed for the high-speed digital computer 
FERDINAND. 


2. The basic equations 
(i) T'resca and von Mises yield conditions 

Referred to cylindrical polar coordinates (r,@,z) the curved bar (Fig. 1) 
is bounded by the circular ares r = a, r = 6; the radii 6 = + and the 
planes z — +h. The applied bending moment M has the effect of bending 
the plate in its own plane. 


Fic. 1. Geometry of the bar. 


Plane-stress conditions are envisaged in which it is assumed that the 
only non-zero stress components are a, and og. These stress components 
are assumed to depend on the radial coordinate r alone so that the equation 
of equilibrium in both elastic and plastic regions is 


do, 1 %—% 


— (), 
dr — r ' (1) 


in the absence of body forces. 
Denoting by u, v the radial and circumferential displacements, respec- 
tively, the strain components in this system of coordinates are 


cu u lev lou cvev 
’ € “7 — —=» + —- - . 

cr Or £20 had roo or r 

together with the longitudinal strain ¢,. Since the shearing stress 7,, is 

zero in all the cases considered here, the shear strain component y, 9 is 


zero everywhere. 
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{ 

The elastic and plastic components of strain will be distinguished by 
1 


a single and a double prime respectively so that 


«-=e¢,+e, ete., (3) 


in a plastic region. Since the stress o, is zero the elastic strain components 
are given by Re ony 
Eeg = og—Vvo, (4) 
ve, v(o,-+-o9) 

where EF is Young’s modulus and v is Poisson’s ratio. 

During the course of the deformation of the plate four possible cases 
are found to occur depending on the applied bending moment and 
distinguished by the number of elastic and plastic regions present in the 
plate. In the first of these the entire plate is elastic, and the components 
of stress and displacement may be obtained from those given in Timo- 
shenko and Goodier (6) for the corresponding plane-strain problem. The 
plate first goes plastic at the inner curved boundary so that in the second 
case there are two distinct regions. In the inner region, a < r < p (Fig. 1), 
the plate is plastic and in the outer region, p < r < b, the plate is elastic. 
This case ends when a second plastic region develops at the outer curved 
boundary. In the third case there are three distinct regions, with the 
regionsa < r < pandé < r < b plastic, and the region p < r < € elastic. 
In the final case the whole plate is plastic and the boundaries r = p and 
r — € coincide. 


In all the four cases the stresses must satisfy the conditions 


o, =0 (r=aandr = Db), 
b 
| rog dr M, 
a 
so that the total applied bending moment is 2Mh. At an elastic-plastic 
or plastic-plastic boundary the normal stress o, must be continuous, 
Material on the boundary of an elastic region must be on the point of 
yielding so that the appropriate yield condition must be satisfied. In view 
of the continuity of o, this condition implies continuity of the circum- 
ferential stress og at an elastic-plastic boundary in all the cases considered 
here. There is a discontinuity in this stress across the plastic-plastic 
boundary of case four. 
The conditions to be satisfied by the displacements are 
u =" _ 0 atO=Oandr 
er 
where k(a-+-b). 
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These conditions imply that the point r = c, @ = 0 remains fixed in space 
and that the radial line segment passing through this point does not rotate. 
Across an elastic-plastic boundary, u, v must be continuous. 


(ii) T'resca yield condition only 
Since all shearing stresses are zero, the principal stresses are o,, og, and 


o,(= 0). The forms of the Tresca yield condition and its associated flow 
rule required here are 


og = —2k, (€", €), €2) = (0, —A,A), (9) 
og—o, = 2k, (€, &, €) = (—A,A, 9), (10) 
where k is the yield stress in shear, dot denotes differentiation with respect 
to time, and A is a non-negative quantity. Equations (9) correspond to 


§% 





2 8X/3—__-o 


| 

















7/6 
Jon 


Fic, 2. The yield conditions. 





the side CD of the yield hexagon (Fig. 2) and equations (10) to the side 
AB of the yield hexagon. In any plastic region governed by equations 
(9) the plastic component of the longitudinal strain is non-zero, so that 
thickening of the plate is to be expected. This thickening is unlikely to 
be important whilst the plastic strain components remain small and is 
neglected here. 


(iii) Von Mises yield condition only 
When o, and og are the only non-zero stresses the von Mises yield condi- 
tion, shown in Fig. 2, has the form 


f = ofa, 04+ 05—3Y* = 0, (1) 
Z 
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where Y is the yield stress in shear. The plastic strain components are 

then given by ef 
A- A(20,— a9) 

d (12) 

of 
A A(2a9g—<a,) 

cO4g 
where again A is a non-negative quantity and dot denotes differentiation 
with respect to time. Equations (12) can be combined to give 

€,(209—<a,) = €4(20,—<59). (13) 

The remaining equation for the plastic strain components is the plastic 
incompressibility condition which has the form 


” 


€ + ej+e, 0. (14) 
This equation also holds for the Tresca yield condition, and is implied 
by equations (9) and (10) in that case. Since small strains alone are con- 
sidered, equation (14) may be integrated to give 


” ” 


ce +egt+e, 0. (15) 


\s for the Tresca yield condition, any tendency of the plate to thicken is 
neglected. 
3. The elastic solution 
Within any elastic region of the bar equations (1), (2), and (4) lead to 
the results 
o, A/r?+ B(L+2logr)+2C 
A jr? + B(34-2logr)-+2C 
$y) B(L--logr)+-C} 
H sin @-+-K cos @4 
A 
E\ r? 
» H cos 80— K sin @+- Fr+4Br0/E y, 


(1-+v)+-2B(1—v)logr— B(1+v)+2C(1 »)| 





where A, B, C, F, H, and K are constants. These equations may also 
be obtained by adapting the results to be found in Timoshenko and 
Goodier (6) for the corresponding plane strain problem. 

For a completely elastic plate equations (16) together with the condi- 
tions (5), (6), and (7) give 
4M ja*h* 
N \ r? 
4M 
N | 


log(b/a) —b? log(b/r) —a? log(r a)| 


24,2 
u*) : , log(b/a)—b? log(b r)—a? log(r/a)!| 
r2 


(b? 


9 
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My, b?— a?) — 2h? log(b/r) — 2a? log(r/a)} 
EN u , : ' 


K, cos 6- 


4Mr{a*b? 
EN\ r? 
+ b?(1—v)log(b/r) + a2(1—v)log(r a), 


(1 -+-v)log(b/a)+- (6? —a?) 4 


oy (b?—a®)rd 


where N — (b?—a?*)?— 4a*b*{log(b/a)}*, (18) 


K, sin 0+- 





and K, is determined by the condition that at the point r — c, 0 = 0 the 
radial displacement is zero. 

For either the Tresca or von Mises yield conditions the plate first 
becomes plastic at the inner curved surface, r — a, with og = —2k or 
a9 Y \3. This was checked numerically for 3/2 < b/a < 100, and is 
almost certainly true for all ba. With b/a — 2-0 the stress distribution 
within a curved bar for Q — 0-129, where 

@ — M)/2ka* — M Yatv3, (19) 
is indicated by one of the curves of Figs. 3 and 4. This value of Q is that 
for which the inner curved boundary is just plastic. 


4. Stress distribution for the Tresca yield condition 
(i) First elastic-plastic stage 

The plate first becomes plastic at the inner curved surface, r — a, when 
ag — —2k and o, — 0. This situation corresponds to the point D of the 
yield hexagon (Fig. 2). It is to be expected that as the bending moment 
M is further increased an inner region of the plate, a < r < p, will become 
plastic. One possibility is that within this region og 2k, o, = 0; but 
the equilibrium equation is then not satisfied. The second possible yield 


condition is 2] 
Oo, — 9 2k, 


so that using the equilibrium equation and one of the conditions (5) we 
have o, 2k log(r/a), o% 2k{1 + log(r/a)}. 

The second of these equations violates the condition (og! < 2k. This 
leaves the final possibility, that the yield condition corresponds to the 
side C'D of the yield hexagon with 


a — 2k, 2k <a, < 0. 

This yield condition together with equation (1) and one of the conditions 
(5) leads to the results 

2k(1—a/r) | 


, 20 
2k j (20) 
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. 3. Radial stress distribution for b/a 2-0 with the 
von Mises yield condition. 
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Circumferential stress distribution for b/a 2-0 with the 


von Mises yield condition. 
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for the non-zero stress components within an inner plastic region, 
a: r < p- 

The derivation of the equations for the stress components within the 
elastic region, p <r < b, using equations (16) together with the appro- 
priate boundary conditions is straightforward and results in 


»]- 42 
| “P[p—a +a log(b/p)|—p(p—a)— 
N,\ r? 

— pa log(r/p)—b?(2—a p)log(b/r)| 
Sye 2 
| 6 Pip —a-}-a log(b p)| 4 2h? — p? - 
N, r? 


09 





— ab? /p—b?(2—a/p)log(b/r) — ap log(r/p)| 
where N, = 2b? log(b/p)—(b?—p*). 
The applied bending moment is found to be given by 


k | 


M . 
N,\ 


(b? — p*)(b? —- a?) — = (b? — p?)? — 2b? log(b/p)[ ap log(b/p) + (p? -aty]. 
-p 
This stage ends when a second plastic region develops at the outer 
curved boundary, r = 6, when og = 2k, o, = 0. 


(ii) The second elastic-plastic stage 
In this stage of the deformation there are three distinct regions of the 
plate. An inner region, a < r < p, satisfies the yield condition (9) so that 
the stress components within this region are given by equations (20). 
Arguments similar to those already used in connexion with the inner 
plastic region indicate that the yield condition appropriate to the outer 
plastic region, € < r < 5, is that of equation (10), corresponding to the 
side AB of the yield hexagon of Fig. 2. Equations (1) and (10), together 
with one of the conditions (5) result in expressions for the stress com- 
ponents in this region of the form 
0, 2k log(b/r) (23) 
og = 2k{1—log(b/r)} 
The central region of the plate, p < r < €, is elastic and equations (16) 
give the expressions 


k 
ga_p2'P(p +a) €2/r?— 1) —(3—a/p)(2—p*) + 2(pa + €*)log(r/p)} 


: jal plp+ay(g? 2—1)—(1 —a/p)(€2—p*) + 2(pa-+ €*)log(r/p)} 
(24) 


Sa 
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for the stress components. The radii of the elastic-plastic boundaries p 
and € are connected by the equation 
2(£* + pa)log(/p) — (€%—p?){3—a/p—2 log(b/€)}, (25) 
and the equation for the bending moment is 
M — ki log(b/p)—(1 —a 2p)(€* — p?) +- 4(b? — &) b__q?)}, 
This stage of the deformation ends when p and € coincide and the whole 
plate is plastic. 
(111) Com pletely plaste stage 
In the final stage of the deformation the whole plate is plastic. This 
case has been considered by Shepherd and Gaydon (5) but the results 
are included here for completeness. Equations (20) hold in the region 
a- r<— p and equations (23) in the region p <r <— 6b. The equation for 
dos l—ap log(b p), 
and the bending moment is given by 


M hkih? +. 2a* — p*— 2pa}j. 


5. Displacements for the Tresca yield condition 

The displacements and longitudinal strain component within an inner 
plastic region are determined by equations (2), (3), (4), (9), and (20) with 
y,9 zero. Substituting from equations (20) into equations (4) results in 


Ke, 2ki(1—v)—a/r} 


Re, 2ki(1—v)+-va/r} | (26) 


Eé, = 2kv(2—a/r) 


for the elastic strain components. The equations for the plastic strain 
components which result from equations (9) are 


a ee 
P 0, age, 0 


r , 


and these equations integrate to give 
0, ene. 0, 
Equations (26) and (27) combine to give 
eer pith 
Ele, + «9+ €,) 2k(1 — 2v)(2 


which together with 
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determine the displacements and the longitudinal strain component. The 
solutions of these equations required here are 


u — H,sin@+ K,cos0+ A,—(2k/E){(1—v)r —alogr} 
v H, cos 6— K, sin 04 Fyr+ Byrd 
A, , 2k{ a\ a | 
3 4 2vi 1 l 
B, ~ + Hl ( } 7 OBr 


Ay (logs ty) 





B, 


where A,, B,, Fy, H,, and A, are constants. These particular solutions are 
chosen so as to have the same general forms in the variable 6 as the 
corresponding functions have in the elastic solution (16). In view of the 
conditions implied in equations (9), the condition €) > 0 must always be 
satisfied within this type of plastic region. 

Using a method similar to that above, equations (2), (3), (4), (10), and 
23) give the equations 


Cu u lov 2k(l—v 
—— f- — — - : 9 Log(b r)—l}, 
or r reo E 

lou cove 

joes 0, 
rover r 
from which the components of displacement within an outer plastic 
region can be determined, The required solutions of these equations are 
D, 2k 


H, sin 64K, cos 6—4 B,r-+ my —F (1—v)rlog(b/r) | 


v — H,cos0— K,sin@+ B,ré4-F,r 


where B,, D),, F,, H,, and K, are constants and the longitudinal strain 
component is given by 


Ee, = Ee, = 2kv{2log(b/r)—1}. (30) 


~ 


The condition to be satisfied in this type of plastic region is (by equations 
(10)) that €) > 0 with 


«9 = 9—«9 = 4B, + D,/r?—2k/E. (31) 


(i) First elastic-plastic stage 

The displacements corresponding to the stresses given by equations 
(20) and (21) for the first elastic-plastic stage are obtained by applying 
conditions (7) to equations (16) in the elastic region, p <r < 6, and 
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equations (28) in the plastic region a <r <p. This operation readily 
gives 
2k| ve 


:. «a | 4k (? Mia 
2p +.— log as 1 }{62(2— —pa' 
E\ ae °B'P “J EN,\r } alla. 


€ 


»)j. 
K, cos 6 = ‘(l—v)r 


E +-va+-a log(p/r)}— 


Ee {b?(2—a p) — pa} 
fj. " 
7 log(p/r)+ EN .. ; 1)o%(2—a pa) 
a {(1—a_ 2p)(b?—p?)— pa log(r/p)—b?(2—a/p)log(b/r)} 
¥a¥y 





2kr (b? 
K, cos @ ee e Pil +-v)[ p—a+a log(b/p)|+ 


aN, r* 


-(1—v)b?(2—a p)log(b/r) + pa(1—v)log(r/p)+ 





+-p(1—v)(p—a)+67(2—a p) -pa| 
(p < <8), 
4k 
” — K,sin 6+ -_. {b>(2—a/p)—pa}ré, 
EN, 
where V, is defined by equation (22) and the constant K, is determined 
by the condition that the radial displacement is zero at the point r = c, 
6 0. 


(ii) The second elastic-plastic stage 


The displacements for this stage of the deformation are obtained by 
a method similar to that used above. Equations (28) and (16) hold in 
the inner plastic region and the elastic region respectively and equations 
(29), (30), and (31) in the outer plastic region. Conditions (7) now give 
4h( ot +- E*)/p 2k a a 
eran (1) + 5 o(2 } +“ log(p/r)! 
E(&—p*) \r E\ \ r r 
. thp( pa + &* 2k 
K, cos 0— —P\P = -{(l—v)r+va+alog(p/r)} 
E(& p*) E 


2k th: + £2) 
~~ log(p r)-4 (pa-r$ ‘ts i 
Er 


E(€*— p?) 


= 
ng eta + €*)log(r/p) —(1 —a/2p)(€?—p?)} 


i a ry i p(p--a)(1- (1) 
w\¢ p”)| r- 


-+-2(pa+-&*)(1—v)log(r/p)—(1—a/p)(€*—p?)(1—v) 4 
+ 2v(£*— p*) + 4(pa-+ €*)| 
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2kv{2 log(b/r)—1} 
a 
E \r? 
atoe-ta(_€) 
E(&— p?) r2 


K, cos 6+- 


—(1—yv)log(6/r) 


ns 
rT '}\ €<r<b), 


| Be) 


sy 4k(E*+ pa) 
- — K,sin 0+ —-._—"— r 8; 
' E(£— p?) 


where A, is determined by the condition u = 0 when r = c, 6 = 0. 


Y -: 


(iii) Completely plastic stage 

Expressions for the displacements within a completely plastic plate may 
be obtained by taking equations (28) in the inner region and equations 
(29), (30), and (31) in the outer region. A continuous solution of the type 
given for the other cases cannot be obtained for this case. The conditions 
implied by equations (9) and (10) cannot be satisfied simultaneously by 
a continuous solution. Since the plate never becomes completely plastic, 
but always contains a thin inner elastic region, a solution for this case is 
to be expected as the limit of the second elastic-plastic case as p — €. 
The first of equations (32) implies that within the inner plastic region the 
longitudinal strain becomes very large as p - € so that some thickening 
of the plate is to be expected. The components of displacement also 
become very large as p-> €. The solution given here is valid for values 
of the displacement small compared with a. Since the magnitude of the 
displacement is largely determined by the factor 


ka 


E(&—p?) 


this solution will be valid for 7-41. With Ek = 400, ba = 2-0, 
pa-— 1-4, &€a = 1-525 it is found that 7 = 0-007 and the value of the 
applied bending moment is then within 99 per cent of the fully plastic 
bending moment. Thus for all practical purposes the solution given here 
is likely to be adequate since some kind of instability will almost certainly 
occur before this value of M is reached. 

It remains to be shown that the conditions to be satisfied by €2 and 
€) in an inner and outer plastic region respectively are fulfilled. Differen- 
tiation with respect to time may be thought of as differentiation with 
respect to p, the radius of the elastic-plastic boundary. 

When the differentiation is carried out the conditions are found to be 
satisfied and the solution presented here is complete. 
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6. Stress distribution for the von Mises yield condition 
(i) The first elastic-plastic stage 

The elastic solution of section 3 ceases to hold when og = —Yv3, 
o, — 0 at the inner curved surface, r = a. A region of the plate near the 
inner boundary becomes plastic, and in this region the two non-zero 
stresses o,, og are determined by the equation of equilibrium (1) and the 
yield condition (11). For the inner plastic region, a < r < p, it is con- 
venient to write the stress components in terms of a parameter ¢ (see 
Fig. 2) defined by 

a, 2Y sin(d— zr), a —2Y sin(d+- 47), (33) 

with 47 <4 < 4a here. The expressions (33) satisfy the yield condition 
(11), and when substituted into the equation of equilibrium (1) result in 


dd cos(¢— bm) cos d/r. 
Ir 


This equation can be integrated to 
r? — a*y* sec dexp(dv3), (34) 
where y == (34/v2)exp(—2v3/12), (35) 
and one of the conditions (5) has been used. The stress components within 
the inner plastic region are determined by equations (33) and (34). 
During this stage of the deformation the region of the plate, p <r < b, 
is elastic. The stress components in this region are readily determined 
from equation (16). Taking ¢ = a on r = p it is found that 
ie log(r/p)— x( 4 v3sin af C08 a " | 
‘ . (36) 
ry log(r p)-+ x(!- 5) V3 sin a— - 5 C08 a | | 
where P, = b?v3sin a—p? cos x 


:, (37) 
p® = a*y* sec aexp(av3) } 
and J, is defined by equation (22). The applied bending moment is given 
by 
| P, 9 ‘ o_: 
M Y LN, (b? —p*)log(b/p) —4v3 6? sin a— 
—p* cos «log(b/p)+ 4v3 a? ha*y?I,(x)}, 

in which /,(«) is defined by 


a 
I(a) = ( tan dexp(dv3) dd. (38) 
in 
The equations given here for the stresses hold until the plate goes plastic 
at the outer curved boundary, r= b. 
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(ii) T'he second elastic-plastic stage 
In this stage of the deformation there are three distinct regions present 
in the bar. An inner region, a < r < p, is plastic with the stress com- 
ponents given by equations (33) and (34). An outer region of the plate, 
& < r < 6, is also plastic and the non-zero stress components in this region 
are determined by the equation of equilibrium (1) and the yield condition 
(11). For the outer plastic region it is convenient to write 
o, = —2Y sin(+- 47), og = 2Y sin({7—y), (39) 
with — ha < & < }r here (see Fig. 2). These expressions satisfy the yield 
condition (11). Equation (1) takes the form 
noah 4 hn) 4 cos yb 
dr r 
which may be integrated to give 
r? — b*y*? sec pexp(—wv3), (40) 


where y is defined by equation (35) and one of the conditions (5) has been 


0, 


used. 

The central region of the bar, p < r < €, is elastic and equations (16) 
determine the stress components in this region. Taking ¢ — «a when r = p 
and 4 = B when r = € it is found that 


2 2 
: logtr/e)—(1 °) + y |? cos x—V38in a 


YP, {, 
&— p2|" ryj |r 
Ps fy 

| 


Cd, 7] 


| 
J 
2 
3) 2 log(r/p) H(1- °)| _— YP cos x+-V3 sin | 


p 
2 
p r 


& 


where » 
P, 


J oor 
- p* cos a+ £2 cos B | 
p? = a*y* sec wexp(av3) 
£2 — h*y* sec Bexp(—Bv3) 
The equation connecting p and € is 
V3(€?—p?)(sina—sin B) = 2P, log(€/p), 
and the applied bending moment is given by 
M = 4Yy*{(a?+-b?)e7*96 —b?],(8)—a?I,(a)}— 

— Yp*{cos « log(é/p) +-4V3(sin «a—sin f)}, 
where J,(«) is defined by equation (38) and J,(8) is given by 

B 

1,8) = | tanpexp(—pv3) dp. (43) 


in 
This stage of the deformation ends when the boundaries r = p and 
€ coincide. 
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(iii) Completely plastic stage 

The case of a completely plastic plate was considered by Shepherd and 
Gaydon (5) but the results are included here for completeness. Equations 
(33) and (34) determine the stress components in the inner plastic region 


and equations (39) and (40) those in the outer region. The equations for 
x and f (defined as in the previous section) are 
4 
x 3 B+ 4m \, (44) 
a* cos Bexp{v3(a+8)} = b? cosa j 
and the applied bending moment in this case is given by 
M 4Y{4v3(a?+-b?) —y*b7I,(B) —y?a7l,(x)} —43 v3 Yp*(sin a—sin B), 


with p given by one of equations (42). 


7. Displacements for the von Mises yield condition 
The components of displacement and the longitudinal strain which 

correspond to the stress components (33) and (34) within an inner plastic 
region are determined by equations (1)—(4), (13), and (15). Substituting 
from equation (33) into equation (4) results in expressions for the elastic 
strain components of the form 

Ke. Y{(1—v)v3sin 6—(1+-v)cos d} 

Re, Y{(1—v)v3sind+(1+v)cos¢} }. (45) 

Ke. 2Vvv3sind 
Since the stress components defined by equations (33) and (34) do not 
depend on the position of the elastic-plastic boundary, equation (13) may 
be integrated to give 

€,(26g—0,) = €9(20,—a5), 
and substituting into this equation from equation (33) we find 
€, cos(d— $7) + € cos(d+ 47) 0. 

Equations (3), (15), (45), and (46) may be combined to give 
y 
€, cos(hd— hz) +-€, cos(d+ dz) FP 
2Y é 
“ 


(1—2v)sin 2¢ | 
(47) 


(l1—2v)v3sind | 


These equations, together with the condition of zero shear strain and 
equations (2), determine the displacements and longitudinal strain in this 
region. 


Solutions of these equations may be obtained by assuming that the 
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circumferential displacement, v, has the same form in this type of plastic 
region as it has in an elastic region, namely, 
v = B,6r+-F,r+ H, cos 0— K, sin 8, (48) 
where B,, F,, H,, and K, are constants. Zero shear strain implies 
Leu, _v 
rob ér r 
so that the radial displacement has the form 
u = H,sin6+ K,cos0+-f(r), (49) 


where f(r) is a function of r alone. The first of equations (47) now gives 
an equation for f of the form 


Feosid— in) +LZeos(d-4 kv) = — B, 008(¢+4+- $7) — pl — 2v)sin 2¢, 
so that, in view of equation (34), 
4 cos ¢ +f cos(d+ 47) 
—ay sect exp(v3, 2) B, cos(¢+ $7) + : (1—2v)sin 24) , 


This equation may be integrated to give the result 


oo r| Det —4 By +4 Bye Hh (4)— (1 —2v)sin(b— 4) |, (50) 


where D, is a constant of integration and J,(¢) is defined by equation (38). 
Substituting for f from equation (50) into equation (49) now gives for 
the radial displacement 


u = H,sin 6+ K, cos 6+ 
4 r| Dye #948, + 4B eK (8)— 5 (1 —20)sin(b— 4) - (51) 

The second of equations (47) gives an expression for the longitudinal strain 
of the form 

€. — 7, (12vy3sing— 

—sin ¢ see(¢—47){ D, e-¢*34+-4B, + 4B, e931, (4)}. (55 
The plastic component of circumferential strain is given by 
<j = €—& 
= De 34-4B,+ 4B, e-%31,(p)+ * cost — $7). 

Equations (33) may be combined to give 


20g—0, = —2YV3cos(¢—47), 
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and since 47 < ¢ < 47 here, this quantity is always negative. The 
quantity A in the second of equations (12) must be a non-negative quantity, 
and reference to equation (54) shows that a sufficient condition for this 
to be the case is that €) < 0 in this region. 
By a method. milar to that used in considering the inner plastic region 
it may be shown that the equations 
€, cos(yh + far) + €g cos(yb— bz) : tl 2v)sin(2yh), 


oYV 


€,-+€9+ € EB (1—2v)V3 sin yf, 
together with the condition of zero shear strain determine the components 
of displacement and the longitudinal strain in an outer plastic region. 
The solutions of these equations required here are 


u H, sin 64 K, cos 64 


N r| Dyexp(y3) 4 B,4 4B, et 31, (pb) mt 2v)sin(ys 4 In) 


4 


H, cos @ K, sin @ } B,r64 Fr 


ml 2v)v3sin bo 





i sinydsee(ys \r){ De?s3 } LB, + A Bieh 1 (b)} | 
where B,, D,, Fy, H,, and K, are constants and J,(y) is defined by equation 
(43). The plastic component of circumferential strain is given by 

. ;, ' Y 

€) D, e844 Bo+-4 Bye! 1b) A V3 costly) +- dar). (56) 
In an outer plastic region the quantity A appearing in the second of 
equations (12) must be non-negative. The expressions (39) may be com 


bined to give ; 
E 209-90, 2Y v3 cos(y+ 47), 


and since this function is positive for the range of values of % considered 
here, A will be non-negative if €) > 0. 

Hill (7) has considered the equations which result when the von Mises 
yield condition and the associated flow rule are assumed in plane-stress 
problems. In terms of the parameters introduced here, thickening or 
thinning is unlikely to be important for any of the values of ys which 
occur and is only likely to be important for ¢ > 42. This value of ¢ occurs 
only for ba > 3-345, so that in practice thickening is unlikely to be 
important, 


(i) First elastic-plastic stage 


The displacements and longitudinal strain for this stage of the deforma- 
tion are determined by applying conditions (7) to equations (48), (51), 
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(52), and (53) in the inner plastic region, a < r < p, and to equations (16) 
in the outer elastic region, p <r < 6. The procedure is straightforward 
and results in 

Y — Y., 

pl -2v)v3 sin d- Boin psec(— hr) x 


aed + = e- #4 1,(¢) —1,(x)|—e* pal +3 cos(a— in)]| 
) | N, N, 1 1 N, | 


Yr{[2P, 


2P, 
E\N, N 


ae 


K, cos 0 + V3 cos(a bn) fom ow3 4 + 


f (a) ~ 1,(¢) je-*3 4 (1 — 2v)sin(¢— in) 
“vy 


Y (2P, 2P, 

pl nen OT Le) — hae 4 

4 4 - 1 

+ V3[cos(a— Ar)e*-%3— cos(d 47)]} | 


r< p), 
4v¥ (P, 


= / a je . | ) 
; iw, 08(r/e) §V3sin « 


rY(2Fi —v)| 1—log(r/p)|+4 


K, cos @ 
: E\N, 





? 2 2 
Py (1 1? ) +(1+-»)?. cos a+(1—v)v3sin a! 
N, r? a 


j 
(p <r <b), 

,’p 

vatd, 

EN, 

where P, and N, are defined by equations (37) and (22) respectively. The 

constant A, is determined by the condition that the radial displacement 

is zero at the point r = c, 0 = 0. 


v K, sin 6+ 


(ii) Second elastic-plastic stage 
The displacements and longitudinal strain in this stage of the deforma- 
tion are again determined by applying conditions (7). In the inner plastic 
region, a < r < p, equations (48), (51), (52), and (53) hold, in the central 
elastic region, p <r < €, equations (16) hold, and in the outer plastic 
region, € < r < b, equations (55) and (56) apply. It is found that 
Y _ Y sin 
€, p" 2v)V3 sin d— a er 4 
< {2P, + 2P, 9 1,(p) —1,(«)]—[2P, + V3(g?— p*)eos(a — gz) je 99} 
ae 2) {[2P, + v3(€2—p?)cos(a— 47) jer-o*3 4 
+ 2PJ1y(x)—1(d)]e-*9-+ 2P, + (1 —2)(€2—p2)sin($— 4} 


K, cos 6 
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ihe Gi (2 2Pfeo -v3__ 1]+ 2P[1,(«)—1(¢)]e -$v3 
ae olf 
+ V3(é*—p*)[e2-*3 cos(a or i —— 


<rc<p), 
4Yp 
— - {P,1 —4V3(€2?—p? 
E@—p)' P, log(r/p) — 4v3(€?—p?)sin a} ) 
Yr 


K, cos 6— ke | 2P,(1—v)[ 1—log(r/p)|+ 


+ P,(1+ Fit 1)+ -(€?—p? rt (149) COs a+ wap: } 
(pxr<é 


, Y sin yb 
a 2v)s <3 
Pa 2v)' iny+ E(é—p?) cos (b+ yn) * 


x {2P,+ 2P, e{ 1,(ys) — 1,(B)]—e-®*3[ 2P, — V3(é?— p?)cos(B+ 47) ]} 
. Y 
K, cos 6— Re _ a ie BW 3[ 2 P, — V3(€2—p?)cos(B+ 47) |— 
. 2 2P, e5{ I,(B)—Iq(y))]-+2P,+(1 —2v)(€2—p2)sin(-+4)} 
= f 2 e-Bv3] 1 OF P, obs T. of 4: 
Rep PAL [Zb)—1,(8)] 
+ ¥3(€2—p?)[e4-8*3 cos(B+ 47) —cos(-+ dn] 
(&<r<b), 





AYP, | P 

E(§?—p*) 

where P, is defined by equation (42). The condition u = 0 at the point 
r = c, 0 = 0 determines the constant K3. 


v = —K,sin 6+ 


(iii) Completely plastic stage 

As with the Tresca yield condition, expressions for the displacements 
which are continuous and of the type given above for the partially plastic 
bar cannot be made to satisfy all the conditions imposed upon them. 
The discussion of the situation given for the Tresca yield condition applies 
here also. 

A numerical check was made to verify that the solutions given in the 
plastic regions of the plate satisfy the conditions on ¢€% for the values of 
b/a considered here (considering differentiation with respect to time as 
differentiation with respect to p), so that the solution given here is 


complete. 
8. Discussion of results 


The integrals J,(¢) and J,() (defined by equations (38) and (43) respec- 
tively) are fundamental to this work, and their values, tabulated at inter- 
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vals of 2° in ¢ and ¥, are given in Table 1 for the range of values needed 
here. 

In what follows a plate of material obeying the Tresca or von Mises 
yield conditions will be called a Tresca or von Mises plate as is appropriate. 
Tables 2 and 3 give the variation of the elastic-plastic boundaries p and 
& with the applied bending moment M for 6/a = 2-0 and 3-0 respectively. 
Values are given for both the Tresca and von Mises yield conditions with 
the same applied bending moment. For both yield conditions the varia- 
tion of p with M in the first elastic-plastic stage is very nearly linear. For 
ba — 2-0 (Table 2) and b/a = 3-0 (Table 3) the maximum difference 
between the values for the Tresca and the von Mises yield conditions 
occurs when the Tresca plate is fully plastic. The difference in p and & 
is of the order of 13 per cent and slightly less than 65 per cent of the 
area of the von Mises plate is plastic when b/a = 2-0. For b/a = 3-0 the 
maximum difference between the Tresca and von Mises values for p and 
€ is of the order of 25 per cent and just under 50 per cent of the von Mises 
plate is plastic when the Tresca plate is fully plastic. These fairly large 
differences in the values of p and & for the two yield conditions may be 
reduced to differences of the order of | per cent for b/a = 2-0 by tabulating 
the variations of p and € with R, where 


R = (M—M,)(M,—™,). 


M, and M, are respectively the values of the applied bending moment M 
when a plastic region first develops and when the plate is fully plastic. 
M, is the same for the two yield conditions but M, is different. This 
variation is given in Table 4 and seems a good way of comparing the 
results for p and € for this particular problem. 

Tables 5, 6, and 7 give the variation of the radial and circumferential 
stresses o, and og with the radial coordinate r for b/a = 2-0 and fixed 
values of the applied bending moment. The agreement of the stresses is 
very close for the three values of M considered here, the graphs of values 
for the two yield conditions being almost indistinguishable. Agreement 
is slightly worse if the stresses are tabulated in terms of R, but is still 
within 5 per cent. The variation of the same stresses assuming the von 
Mises yield condition, with b/a = 2-0, is shown in Figs. 3 and 4 for several 
values of the applied bending moment. 

Most of this numerical work was programmed for the digital computer 
FERDINAND. 

In conclusion, it is seen that in this particular problem the stress distri- 
bution is relatively independent of the yield condition. The position of 
the elastic-plastic boundaries is much more dependent upon the yield 

5092.51 Aa 
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condition, but use of Table 4 may be of use in overcoming this difference. 
It is, of course, difficult to estimate how general these results may be with 
regard to other problems. 
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TABLE | 
Values of the integrals 


¢ ¥ 
I,(¢) | tangexp(dv3)dd and IY) = | tanvexp(—yv3) dy 


iz —itr 


for a range of values of ¢ and 
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TABLE 2 


Values of pa and & a for a range of values of 
M Ya*v3 for ba — 2-0 


Y M 2ka* 





Von A Tresca 





pa éa 











TABLE 


a range of values of 
3-0 


Values of pa and &a for 


( M 2ka® — M Ya*v3 for b/a 





Tresca 





fa 
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TABLE 4 


Values of pia and €/a for a range of values of R = (M—M,)/(M,—Mo) 
for b/a = 2-0. My, and M, are the initial and fully plastic bending 


moments respectively 





Von Mise Tresca 


fa fa 








°° 
“gob 
‘RIG 
738 
O04 
sgh 
“S4 


“458 

















TABLE 5 


Values of the stresses a,, ag for a range of values of ra with 
ba — 20and @ = 0-1598 





Von Vises Tresca 





a, YN3 Oy YV3 





°o 1°o C 
ooho4 orosgs 
o7ogos O0%47 
01439 
o1676 
or1706 
o'1s9g2 
0°1 377 
o*109g2 
o°0759 
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: TABLE 6 
Values of the stresses o,, og for a range of values of r/a with 
ba = 20and Q = 0-1926 





Von Mises Tresca 





Yv3 Og YN3 2k 





"O595 
O1i7s 
1629 
"rg9gt 
2073 
*1Q6o0 
“1701 
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"0953 
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TABLE 7 


Values of the stresses o,, og for a range of values of r/a with 
ba — 20and Q = 02375 





Tresca 
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0°2 308 
0°2857 
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SUMMARY 

The formal integration of partial differential equations defined over regions with 
elliptical boundaries may, in suitable circumstu :es, be achieved by expansion in 
terms of Mathieu functions. The numerical integration of such equations by this 
method, however, would be in general very difficult since numerical information 
in respect of the Mathieu functions, and more particularly their zeros, is extremely 
limited. 

The method described in this paper is intended as a practical alternative technique 
which, in addition, should be particularly suited to the solution of non-separable 
eigenvalue problems. To demonstrate it, investigations of the modes of vibration 
of stretched membranes with constant and variable densities respectively have been 
undertaken. In either case the problem is finally reduced to the solution of a simple 
set of simultaneous algebraic equations requiring relatively little labour to solve, 


1. Introduction 
AN investigation of the modes of vibration of stretched, homogeneous, 
loss-free membranes with elliptical boundaries is well known to depend 


upon the solution, for the eigenvibrations 6 and associated eigenfrequencies 
A, of the partial differential equation 


“2p et9 
ae + a + 4Ah*(cosh 2 —cos 2) = 0 (1.1) 
subject to the boundary condition 
OE5.9) = 0 for O <n < 2m. (1.2) 
Here / is a measure of the distance between the foci of the ellipse and 
£ and » are elliptical coordinates lying in the ranges 0 < € < &, and 
0 < » < 2m respectively. Separation of equation (1.1), putting 
AE, n) = BE»), 
leads to the ordinary differential equations 
d*4 
—(a—2q cosh 2 = 
d*p 
dy? 


(Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 3, 1960) 


and + (a—2q cos 2n)b = 0, 
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where a is a separation constant and 2q = }3h?A. These are Mathieu's 
equations (1), the appropriate solutions of which are the Mathieu functions 
of the first kind of integral order, namely, ce,,(7,q) and se,, ,,(7,q) of (1.4) 
and the modified functions 


Ce,,(€.9) . ce, (t€, q) 
and Se, 1(€.9) _—™ ise, .(t&,q) 


of (1.3). The integral suffix m refers to the separation constant a,, for 
m = 0, I, 3...., whilst ¢ denotes .(—1). Further, the boundary condition 
(1.2) requires that 
( 'Cm(Eos q) 0, 
S€m+i(€o> 9) 0, 


giving, for each m, an infinity of zeros q,,,. Gmsip Where r 
The appropriate solution of the membrane problem is 


AE, ”) > A,, Ce, (€,q) ce, (7,9) T > B,, Sen irlE, V8Cm1(7 9): (1.5) 
) 


m 0 m ( 
where each individual term of the series corresponds to a normal mode of 
vibration when the appropriate value of g has been determined. 

The direct evaluation of the parametric zeros q,,,. Fn «1, for a prescribed 
shape requires a considerable amount of labour so that the practical utility 
of the above analysis is greatly curtailed. A limited amount of numerical 
information in relation to the Mathieu functions is available (2), and in 
a recent paper Daymond (3) has used this to caleulate the eigenfrequency 
of the principal mode for constant density membranes of differing eccen- 
tricity. However such numerical information as there is would be in- 
sufficient to allow a comprehensive investigation of the modes of (1.1) to 
be undertaken without prohibitive labour, and this would be even more 
forcibly the case if. for example, the modes of a variable density mem- 
brane were to be investigated since the governing equation in this instance 
would be, in general, non-separable. 


In this paper we suggest an alternative and more direct technique which 


avoids separation of variables and thus the determination of the separation 


constant which a conventional analysis of the homogeneous membrane 
would introduce. 

Investigations of the modes of vibration of stretched elliptic membranes 
with constant and variable densities respectively have been undertaken 
with a view to illustrating the facility of the method. In the former case 
a further simplification can be made and is suggested as a distinct alterna- 
tive to the normal method of evaluating the Mathieu functions themselves. 
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2. Formulation of the solution 
For an elliptic membrane of variable density p(£,7) the eigenvalue 
equation is 
29 ot 
; ee “ + SAR? p(E, 7){cosh 2 —cos 2n}6 = 0. 2.1) 
oF * ig” 


We postulate a fictitious range —& < € < &, for the coordinate € and put 


g — (= -1) 


y= y 
Using the transformations (2.2) and putting R = 2£,/7 the equation (2.1) 
becomes 
c*0 ,078 Apes , : ak aie 

+ R?— + JAR? RR? F(x, y){cosh R(2x—7)—cos 2y}8@ = 0 (2.3) 
cx? cy 
defined over the rectangular region 0 <4 < 7, 0 < y < 2m, the trans- 
formed function p(€,) being represented here by F(x, y). 

Further, in terms of the new coordinates x and y, we rewrite the original 
boundary condition (1.2) as 

A(0,y) = 0 = On,y) for O<y < 2n. 2.4) 

In general the problem posed by the modified equation (2.3) and the 
boundary conditions (2.4) is readily amenable to solution by a formal 
Fourier expansion in terms of products of appropriate trigonometrical 
functions. Moreover, for our solution to be a proper physical one, two 
further conditions require to be satisfied, namely, 


(i) (47, y) = O(4n, 27—y), 


(ii) lim be {O(x,y)} = — lim “O(a, 2a y)}, 
r x 


rin ©. rin 


ensuring continuity of displacement and of the gradient, respectively, in 
crossing the interfocal line orthogonally. 
Initially we propose as a solution an expansion of the general form 


° cos n 
Ox. y) = > Sa,,,,sinmex . y 
m on ; sin ny 


2.5) 
which implicitly satisfies the boundary condition (2.4). The choice of 
either factor cosny or sinny in the expansion (2.5) will be governed in 
any instance by the geometrical features of the mode under consideration. 
Also the appropriate values of the integral suffixes m and n must be 
assigned in each case. 

Each basic mode may be classified in one of four categories according 
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as it exhibits symmetry or antisymmetry about each of the axes of the 
ellipse. Thus the basic modes are each characterized by one of the follow- 
ing characteristics: 

(a) symmetry about both axes; 

(6) antisymmetry about the major axis and symmetry about the minor 
axis; 

(c) antisymmetry about the minor axis and symmetry about the major 
axis; 

(d) antisymmetry about both axes. 
Expansions satisfying the properties (a), (6), (c), and (d) respectively as 
well as the continuity conditions stated above are readily found to be 


= ¥ ¥a,,., sin mz cos ny (2.6) 


mn 


for m a % sn = 0, 2, 4....:; 


> 3 g,,.., sin ma sin ny 


m n 


for m 2, 4,6 oe Par oe See 


> > a,,.,, sin mz cos ny 


m,n 


for m 


rye 


mn 


sin mx sin ny 


m,n 


fc rm 


3. The homogeneous membrane 

Substituting the expansion (2.6) into the equation (2.3) with F(z,y) = 1 
and writing 1h2R2A — A (3 1) 
we have immediately 


> ¥ (m?+ R*n?®)a 


mn 


sin mz cos ny 


mn 


A{cosh R(2x—7)—cos 2y} ¥ }a,,,, sinmaxcosny. (3.2) 
mn 


We may now expand cosh R(2x—7) as a half-range Fourier cosine series 
in 2, i.e. 


t 2 
cosh R(2r—z7) — 94. = > 5, cos qx, (3.3) 


7 7 @ 


4Rsinh Rx , 
where b= 4k (q = 0, 2, 4, 6,...). 
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Also the familiar orthogonality properties of the trigonometrical func- 
tions imply that for gq + 0 


(q = |p—m)) 


jn g=p-—m 
(3.4) 


: ' jr qg=m—p 

| cos gx sin px sin mx dx = 
: —}7 q=m+p 
0 otherwise. 


Therefore, substituting (3.3) into (3.2), multiplying across the resulting 

equation by sinmzcosny for each m and n, and integrating over the 

rectangle 0 <2 <7, 0 <y < 2m and using (3.4), we obtain a doubly 

infinite set of linear, homogeneous, algebraic equations in the elements 

a,,, and containing the parameter A. A typical equation of the system is 

(1, ,—-AJ™*)a,,,, = A| > Ima Ka 
v 


mn m,n pn “pn 


m,n q— Tay,» +3) (p+ m), (3.5) 
wae = m*+ R*n?, 


1 (b 
7 


p-m —b, 1m): 


and in general K = L =}. 


However, when n = 2 the value of the factor K becomes unity, Z remain- 
ing at }. 

The system of equations typified by (3.5) can be readily solved by 
numerical methods. The system is homogeneous so that it provides non- 
trivial solutions for the elements a,,,, only for certain values of the para- 
meter A. That is, each solution gives an eigenvalue A together with the 
coefficients of the corresponding eigenfunction of (2.3), thus providing the 
required information relating to those modes whose geometry satisfies 
the expansion (2.6). 

Following the above procedure sets of equations similar to the set 
typified by (3.5) are deduced when the expansions (2.7), (2.8), and (2.9), 
respectively, are substituted into (2.3). For example, when (2.9) is em- 
ployed, the only difference in the form of the equations obtained is that 
the value of the factor K is universally 4. Thus we shall refer unam- 
biguously to the equations (3.5) since they are typical of the equations 
relating to each type of mode. 


4. The variable density membrane 
We consider in this section a density function p(€, 7) given by 


p(é.n) = ( _cosh®é cos?y _ sinh®é int) 


cosh,  sinh*, ) 
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The eigenvalue equation (2.1) incorporating this expression for p(£, ) is 
one which arises in laminar forced convection when a heated fluid flows 
through an elliptic duct with an established velocity profile which is 
proportional to p(é, 7) as defined in (4.1). The boundary condition (1.2) 
further requires that the walls of the duct be maintained at a constant 
temperature different from that of the fluid at entry. 
Using the transformations (2.2) and rewriting (4.1) as 
‘osh?2 1 cost nht Rie — ivan? 
F(x.) (1 cosh? R(x see y «sinh R(a -$e)sin uv) (4.2) 
coshé, sinh?é, 

substitution of (4.2) into the transformed governing equation (2.3) gives 
on simplification 


AG(a, y)0 7 (4.3) 


where h? R?X/2 
and 

G(x, y) leosh P( 2x — m){(1—a)-+ acos 4y}— 

cosh 2 R(2x—7){B 4 «cos 2y}— (1 —a)cos 2y + 8 cos 4y], 

with , | /2 sinh? 2£,; B = cosh 2,/2 sinh? 2€,. 

We require a Fourier series development of the above expression for 
G(x, y) and in this connexion we use (3.3) for cosh R(2x--7) and put 

‘ 


» 
cosh 2 R( 2x2 —7) fi a > C, COB GX, 
7 7 @q 


where SK eink 38a (q = 0, 2, 4, 6.,...). 
q+ 16OR? 

From this stage onwards the analysis follows along the lines set out 
above for the homogeneous case. For example, if the expansion (2.6) is 
substituted into (4.3). we eventually obtain a set of algebraic equations 
each similar in character to (3.5). Of this set the following is a typical 
equation (for m 1, 3, 5, 7,... and n 0, 6, 8, 10,...), 


\ 


7 


(Ain Ent Le (Ay in-a + 4pnse) t+ Mpa (a 


’ 
pn? T An ns2ds> 


(4.4) 
where m2 Rn? 


(1 x)(6 pom a Bic pom Cy +m) 


l iW: . 
2 (5, m b, m) T Ba Om,p 


l ’ . } ) 
2 v{¢ pom ¢ p +m) 5( | x) Omp? 





THE FREE VIBRATIONS OF MEMBRANES 365 


and 6,,,, is the Kronecker delta. It will be noted that, in detailing the 
integral values which m and n may assume in this equation, the two values 
n = 2 and 4 have been excluded. The reason for this is that the equations 
displaying a term a,,,, explicitly on the left-hand side for m = 2 and n = 4 
respectively have slight irregularities in their forms compared with the 
typical equation (4.4). Written separately these equations are 


7 ae Ay 


Nye 


2 { +2 | 270, i 
> (Amza, a, ? LM3 an6 t M3 (24,6 7 ay a)i 


m,2 


and 


en 

Am | | 4 ' 

Ina ma > (Aga, a, + Lint (20 hy, o7 4,4) t Mis (4,2 1 1 y,6)§> 
P 


where, with the exception of 


Kms (1 ha)(b — Bism) Ble — Corum) \Br ae 
the factors J, A, L, and M are as defined above. 
Similar sets of equations to the one typified by (4.4) would of course 
be obtained if expansions (2.7), (2.8), or (2.9) were used for modes dis 
playing antisymmetrical features. 


5. Simplification of the method for F(x, y) — | 

It has already been pointed out that each basic mode of the separable 
equation (1.1) is represented by an individual term of the expansion (1.5). 
Thus, in elliptical coordinates, the principal mode is obtained from the 
single product term Ce,(£,q)ceg(y,q). On the other hand, the solution is 
given by (2.6) of the present paper as 


Ax,y) = > da, sin ma cos ny 
mn 
for m = 1, 3, 5,..., mn = 0, 2, 4, 6,..... Therefore if we write 


yee 


Ceg(x,q) = > Ajsinix (i 


ceg(¥,q) = > Bycosjy (j 
} 


then clearly a A,, B 


mon m "” 
and instead of the set of algebraic equations typified by (3.5) we would 
have an equivalent set typified by 


( mn), B, = ALS Im, B,—KA,, B, g—LAy By .2| 
yp 


m,n mn pn 
(p xm), (5.1) 
where the factors J,,,,, J™", K, and L are those defined for the equation 
(3.5), 
This modified procedure would, in general, lead to labour-saving in 
numerical work, particularly if higher modes or comparatively low modes 
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for an elliptical shape with eccentricity near unity were considered. For, 
comparing Table 1 with Table 4 (see section 6) it is easily seen that as 
against some 4mn coefficients to be evaluated in the first case, only (m-+-n) 
need be evaluated in the modified method. Also the table of Mathieu 
functions (2) would suggest that for the modes mentioned the behaviour 
of the coefficients B, would be of the form B, = 1, B, ~ —1, By ~ +1, 
and so on. Multiplying these by the appropriate coefficients A,, it might 
be expected that a considerable number of coefficients a,,,, would need 
to be calculated before a comparable degree of convergence in these would 
be obtained. In such cases, therefore, the saving in labour if the above 
modified procedure were used would probably be considerable. 

The essential merit of this simplification is that in the form of the 
equations (5.1) we have all the advantages of separation without the 
complications that a separation constant would introduce. For this reason 
the modified method described above might usefully be considered as an 


alternative to normal computing practice for evaluating the Mathieu 
functions themselves. 


6. Numerical solutions 


Several functions characteristic of the numerical problem likely to be 
encountered have been selected for computation. First of all, in relation 
to a homogeneous membrane of eccentricity 0-5, we have computed the 
principal mode, i.e. the solution corresponding to the least value of A satis- 
fying the equations (3.5) associated with the expansion (2.6). Secondly, 
these equations (3.5) have been solved in relation to the expansion (2.9) 
to give the first mode characterized by antisymmetry about both axes of 
the ellipse. The numerical results for these functions, including the eigen- 
function coefficients a,,,, and the pulsatance vA, are given in Tables 1 and 
2 respectively. 

The principal mode of the variable density membrane considered in 
section 4 has also been evaluated for an ellipse having major and minor 
axes in the ratio 2:1. The numerical detail for this function is included 
in Table 3. Finally, Table 4 shows the numerical result for the evaluation 
of the principal mode of the homogeneous case obtained by the simplified 
method of section 5. 

The area of the ellipse has been taken to be m in each case so that 


h® = 2 cosech 2&4, whence A is obtained from the computed value of A on 
using the relationship : 
E \ sinh 2€, 


R2 


A 


where the eccentricity is sech &,. 
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TABLE 1 
Values of ay,» 





° 


2 





1°0000 
0°0553 
0°0065 
0°0013 
0°0004 
o"ooo! 
o°ooo! 

















Value of vA = 


2-4166 (2-4166 (Daymond)). 


9 


- 


TABLE 
Values of a 


mn 





4 




















vA = 5-158. 


TABLE 3 


Values of a 


m,n 





4 











—o2189 
00094 
00008 
00002 














vA = 2-983. 


TABLE 4 





Am 








1°0000 
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The solution of the algebraic equations, such as those typified by the 


general equation (4.4), is a standard computation and is attainable, for 
example, by relaxation methods. However, in virtue of the presence of 
the factor (m?+ R*n?), the term exclusively on the left-hand side of any 
equation is entirely dominant. In view of this and the additional con- 
sideration that the eigenfunction coefficients a,, , converged quite rapidly 


mn 
it was found convenient to use a simple method of iteration which has 
been described in (4). 
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SUMMARY 


In this paper the complex variable method has been applied to the solution of 
problems of dynamic plane elasticity for anisotropic media. The stress and displace- 
ment components are obtained in terms of two analytic functions of two different 
complex variables and the solutions are given for a moving parabolic punch and 
a moving dislocation. 


1. Introduction 

ComPLEX variable methods have been applied to the solution of problems 
of dynamic plane elasticity for isotropic media by Radok (1). He obtained 
the stress and displacement components in terms of two analytic functions 
of two different complex variables and applied the results to obtain the 
solutions for (i) a moving Griffith crack, (ii) a moving parabolic punch, 
and (iii) a moving dislocation. Sneddon (2) has used a different method 
for solving such problems. In this paper Radok’s method is modified for 
two-dimensional anisotropic media and applied to the problems corre- 
sponding to (ii) and (iii). 


2. Deduction of stress and displacement components 
It is assumed here that the stress and displacement components o,, 0 
T,y, U, v are connected by the equations 


C cv 
j.¢ 
1 12- 
p cy 


ev | 


vy? 


the equations of motion being 


a ‘ n9 ; ‘ 
Ce, 4 ey - P. u CTry | O 
Cx cy 


. as uv — 
ce cr cy 
where p is the density of the medium. 
Eliminating 7,,, between the equations (3), and using the values of éu/éz, 


(Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 3, 1960 
5092.51 Bb 
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év/éy obtained by solving (1), it is seen that 


- 29 . a2 2 c2 

eo” __ Plates) € - |; _ Plieten) ¢ Je 
ee} ag , Ae) “a aga) u" 
CX* C4 Cog — Cy QC gy Ct CY” Cy, Cog — Cy Cay Ct 

Hence o,, 0, may be written in the form 


o, = E = ple 4 C12) « v| 
cy* ot 


g2 
31 C22 — (12 21 ; 
— oe ; (Cor + C2) “le ) 
. OX* Cy, Cog—C yg yy OF? 
where U is a function of x, y, t. Again the equation of compatibility may 
be written in the form 


(4) 


2 a2 he 
Sus Cos) . (2eesGCead : p(c Cc . o 
(C41 Cog C12 Coa — = 66 29) Hg T (=O22 Ces) =; 22 a1) 5 Tt 
_s ” 2)" ex? cy? ot? sd 


2 2 2 
oo € € 
De . ~ . » — & 
(2 1166) » at (11 Cop — Cre Ca “60 12) oa p(Cu- C12) a8 ao, = 9%. (5) 


It is assumed that x and ¢ occur only in the form x«—ct, where c is a 
constant. Hence if we set = xr—ct, » = y, we have 
9 “9 9 “9 a) 
co c* 1 oc Ce a 
x2 682? ft?’ cy? én?’ 
and using these results and substituting for ¢,, o, from equations (4) in 
equation (5) we find after some simplification that 
al au au 
{ B ( = 0 (6) 


ef! e&n? Ont 
where re. 
22 “66: 
B = (1, C22—€19Co1 — C66 C12 — C66 C21) — PC*(Co2 + Coe), 
A (pe* C1,)(pe?- Cee): 
It will be assumed that the roots of the equation 
Czr?+ Br+A = 0 (7) 
are —.*, —f* (a, 8 being positive and real), which assumption is valid 
provided that 


(i) P 11 C22 — C12 Car — Coe (Cr2 + C2) 
22 7 Cee 
ii) pe? must not lie between c,, and c,g,, and 
1 66 
(iii) the roots of the equation 
c 9 
x*(Coa—Cog)” 
») > - . . - » 
— 2[ (Coo+ Coe )(C11 C22 — C12 C21 — Coe C12 — Cee C21) — 2022 Coe (C11 + Cee) |+ 
2 “i 
+- [C43 Coo —Cy2 Coy — Coe (Cr2 + C21) |? — 411 Coe Se = O 


must either be complex, or, if they are real, pc? must not lie between them. 
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In the case of wood, using the values for spruce and oak given by Green 
and Taylor (3), it is found that all these conditions are satisfied if c* < Cg¢/p. 
If the above conditions are satisfied, equation (6) may be written 


(steel eget ae)” ~ 


so that (6) is satisfied by 
U = R(z,)+A(z)+ AG) +e), 
where z, = €+ian, z, = £4 1B». 


Writing 
4 “oe . dF, y dis 


as > 
dz, dz, dz, dz, 


one obtains easily 


pcr(ey, +Cy2) \yy 4 fg rts pcr(Cy, + C49) \y’ (9) 


C11 Coa “tratal” \ C11 22 C19 Cay) 


2/ A. 
— PMT [rete +¥ (10) 
C11 ©2212 C21 
ooyT aay? 
since Oe nT on SA 


o€? c2 ot? 


and OU _ _2rela%$’ +BY’). 
én 
Using the expressions for éu éx, év/éy in terms of o,, 0, and integrating, 
one obtains 
» 
—— re[ (¢gg. 0? + Cy + pc?) + (Cg B® + C2 +-pe®)p], (11) 
C1122 — C1221 
» 


~ . , On a (‘ , Cn" 3} 9 
——_— retlic + ¢4 -— — 9 12 
©1122 —Cy2 Can ( ail x B+ B ? ated 
» 
/ 3) ; 2 Cy —pe*| ,, 
ee BOE IC +-(C,9-+Ce, + er pet SE, 
C1, Cag — C19 Cay (| - 1a Cat pe)s om i? ” 


+- {con (Cya + Co, + pe*)B + cu ely. (13) 


The constants of integration have been omitted as they refer only to rigid 
body motions. 


3. Applications 


The above results will now be applied to the problems of (i) a moving 
parabolic punch, and (ii) a moving dislocation. 
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(i) A parabolic punch moving over a half-plane 


The boundary conditions to be satisfied are 


t,, ~ 0, for all points on the -axis, 


a, — 0, for \€| > lon the -axis, 


vy = f(é) = d—&/2R, for \€| < lon the €-axis, 
R being large. 
The solution can be obtained in a manner analogous to that of Radok 
and is given by 
(z,) Ay(2,), ib(Z,) By(2.), 


where 


and 


9 


pe* 


B 


c 


(ya + Ca, + pe?) — (a? + 


12), aa x? + a8 + B)+ (Chat Ce 


(C41 Cog —Cy2C)| » 


Caq(a? + B+ B®) + (Cy2+ a1 + pe) 


where 


+B CoolC11 pe*) ~~ 2) p2 
Bt B («+ B)(a® + B?) 


a " os 
Mo = C9221 9B(x+ B) + (Cy + pe? )(ey, — pe?) 


and «, B are given by (7) 





PROBLEMS OF DYNAMIC PLANE ELASTICITY 


(ii) Moving dislocation 
In this case the boundary conditions are 
o, = 0 | 


for - 0 
an o| ’ 


where [uw], denotes the change undergone by u for one circuit of L. 
Here ¢ and #% have the same values as those given by Radok, namely, 


Ci) =b log 21> yb —hb log 29. 


Hence, 


» 
u . - re b| (Cgy x? 4 C49 + pe log 2, — (Cgg B® +-¢y2 +- pe* log z.], 


2 
—~ Jlog “| 


C11 ©oa~ C12 C21 


” . 2 
r c 
: - wi 11 pe e uN 
t re ib (ca x + : Jlow (rab } 


C11 a2 12 C21 


B 


The results obtained agree with those given by Radok if 
= » 
Cy, = Cag = A+ 2p, Chg = Cy, = A, Cos = B- 
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SUMMARY 
The use of an automatic computer for the direct numerical solution of a cylindrical 
heat-conduction problem is considered. It is shown that accurate solutions can be 
obtained easily and rapidly by use of the Crank—Nicolson implicit method. Attention 
is given to the adequacy of the finite-difference representation in the neighbourhood 
of a singularity at the boundary. 


1. Introduction 

In arecent paper (1) Andrews has given an approximate analytical solution 
of the heat-flow equation for the cooling of a spinning thread-line. He 
remarks that the apparent success of the methods he outlines is fortunate 
in view of the excessive labour of carrying out a direct numerical solution 
of the equation owing to the necessity of using very small intervals in the 
integration steps. 

One purpose of this paper is to point out that the maximum time-interval 
permissible in a numerical solution of the heat-flow equation can be con- 
siderably increased by the use of implicit rather than explicit methods of 
solution. Although the labour in computing a solution of Andrews’s 
problem by implicit methods is still heavy on desk machines, accurate 
solutions can be obtained easily and very rapidly with the aid of an auto- 
matic computer. 

The other purpose of this paper is to discuss the adequacy of finite- 
difference representations of the heat-flow equation in the vicinity of a 
singularity in the boundary conditions. 


2. Statement of the problem 


The equationt solved by Andrews is 


6 2 a 
* : ce l cb (0 4 ,< 1) 
ct Or-+ Tor 


Z 


with the initial condition 


Q 1 forallratt=0 


+ The notation adopted in this paper differs slightly from that used in (1). 


[Quart. Journ. Mech. and Applied Math., Vol. XIII, Pt. 3, 1960) 
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and the boundary conditions 


e = 0 (r = 0), reed = — F(t) (r a 1). (3) 


or 


In these equations, Z is a constant and F(t) an empirical function of t. 
The occurrence of a time-dependent function in the boundary conditions 
(3) prevents the solution of equation (1) from being separable in the vari- 
ables r and t. This makes the analytic treatment more difficult but it adds 
little complication to the numerical process outlined in the next section. 


3. Finite-difference approximations 

We apply the Crank—Nicolson implicit method (2), which is known to be 
stablet for all values of the intervals used in ¢t and r. Let us denote by 
6,, ®;, and ¢; the values of 0(t), 0(t+-48t), and 6(t-+-8t) respectively at the 
points r; = ¢ 6r, where i = 0, 1, 2,..., n and dr = 1/n. Then a finite- 
difference representation of equation (1) is 

se 0 Sanya beer 26+ $; 1 fetes +6; 4 —26;+ 

$6,141 Et _C14,) C09) 

in which C,(¢;), C,(0;), and C(®,) are ‘difference-correction’ terms which 
can be neglected in computing a first approximation. C;,(¢;) and C,(@;) are 
expressible in terms of differences in the r-direction whilst C(®,) involves 
differences in the ¢t-direction. In standard central-difference notation, 
defined, for example, on p. 36 of (3), the leading terms in the expressions 
for C(8;) and C(®,) are 


C(0;) = 


Z 
+ 5p Oo) (4) 


04... 


i? 12 r 


+ 886 
62 


CD) = — 87D; + 997 O;+-..., 


where the suffix r or ¢ denotes a differencing operation in the r- or t-direction 
respectively. 

The term containing C;, in (4) will be neglected throughout the subsequent 
work. Its retention would add considerably to the complexity and storage 
requirements of the process on an automatic computer and it is usually 
preferable to make its contribution negligible by reducing the interval 6. 
This procedure, of course, increases the number of time steps in the calcula- 
tion but does not affect the storage requirements of the programme. The 
situation with regard to the terms C,(¢,) and C;(;) is somewhat different. 
Firstly, their retention is easier in that they are computed from values 


(5) 


+ For a discussion of the stability and convergence of the method see section 6. 
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which must be stored currently inside the machine. Secondly, the require- 
ments of the machine are likely to be more restrictive on the permissible 
values of dr than of dt. The C, term will be ignored in the rest of the paper. 

Let us suppose that at time ¢ the values of 65, 0,,..., 6,, are known. Equa- 
tion (4) with i = 1, 2...., n—l1 then gives a simultaneous set of n—1 
equations for the n+ 1 unknowns dp, 4),..., ¢,,. The two additional equations 
required may be obtained by use of (4) with « = 0 and n, which introduces 
two further unknowns ¢_,, ¢,,.,, and use of the finite-difference representa- 
tion of the boundary conditions (3), which enables us to eliminate ¢_, 
and ¢,.,. Thus at i = 0 is represented by the 
symmetry condition 


a 


0 the condition ¢@/ér 


O_; O44, $..= $4 
whilst the term (1/r)é@ ér is replaced by its limiting value €70/ér?. The 
condition (1 @)é@ ér —- F(t) ati 


(6) 


n is represented by 


6 —§ PS 
nl n-1 —_ p36 | 
26r bor 


¢,, .,—4, ] - 
_ a A —— pd5g,, +... 
zor bor 


—6, F(t) 


Equations (4)-(7) may be combined and rewritten in the form 
Pi; 1 9: %;4 (2—p,)d; 1 
+-1,6;—(2- 


i i 


P; 0, - Pi); at C'(d;) + C(8;) (i — 0, :. - en n), 
where (Sr)? 
Go = 2+2-—- 


\ 


| 
21 
| 


n n 


Equations (8) are a simultaneous set of linear equations to be solved at 
In many instances sufficient accuracy 
will be obtained by neglecting the correction terms C,(¢;) and C,(4;) 


each time step for the unknowns ¢,. 


Co(9%) = 139°05+-... 


aly ZO. 914 45r5rF(t-+ at) 


St | 


wAwsn 


ae 2(1+-36r)dr F(t) 


53, 154 
. hud 6, ; {3° 6, “T eee 
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altogether; if the difference-correction terms cannot be neglected the 
equations may be solved iteratively. Equations (8) are first solved with 
C(¢,) replaced by zero and the correction term evaluated from the approxi- 
mate values of ¢; so obtained. The correction is then inserted on the right 
of (8) and the equations re-solved. The cycle is repeated if necessary until 
consistency is achieved. 

It is worth noting that if the original problem were non-linear the equa- 
tions corresponding to (8) would be non-linear algebraic equations, which 
would usually be solved iteratively. The incorporation of a difference 
correction would then be proportionately less time-consuming. 


4. Numerical solution of the finite-difference equations 
We have to consider 

Pibisi—Ubit(2—P) bin = 5; (t= 0,1,...,n), (12) 
where p; and q; are given by (9), (10), and (11) and the b; are known. Equa- 
tions (12) are a set of (n+ 1) simultaneous equations for which the matrix 
of coefficients is of triple-diagonal form. A method of solving such a set 
which takes advantage of the large number of zero coefficients in the 
matrix is the following well-known recursive procedure. Setting 


X% = do» So = by (13) 


we first compute the sequences 


a ea. ae ae (14) 


= A Reig. 


Then 


1 , 
and d; = —(Pibin—S) (6 = n—1,n—2,..., 0). 
xy 


This technique is equivalent to the familiar process of Gaussian elimination, 
described in (3, p. 4), adapted to take advantage of the special form of the 
equations (12). The multipliers VM, at each stage are given by 
9 
M,=“—"* (§ =1,2....,n). (18) 
yy 
It is not difficult to show that in the present problem M, < 2 and the 
remaining M, < 1, so that no difficulties due to large multipliers can arise. 
It should be noted that for the problem under discussion ag, a,..., &—1 
are independent of the time-variable and can be calculated once and for all 
at the start of the computation. 
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5. Example 


Andrews’s problem has been programmed and solved on ACE at the 
National Physical Laboratory. The required values of F(t) were computed 
from a polynomial approximation to this empirical function. Table 1 
gives the values of F(t) and the solution (r,t) at the points r = 0, 4, and 1 
computed for the cxse Z = 25. The values 6*(r,t) were obtained using 
Andrews’s approximate formula and are given for comparison. 


TABLE 1 





4, t) : O*(4, t) O*(1, t) 





000 “OO . *000 1°000 
“R16 


2812 “712 


712 


I 
o-g18 o "9: “gt 
° 


*HO2 ‘Sol 


° 

° 
"703 0°605 

° 

° 


500 “410 


"332 


266 


410 
"333 
267 212 
“168 
“133 


05 


O° 213 
169 
Ol 34 
o'106 07053 


°90000 




















The values of @(r,t) were obtained by taking intervals of ;, and } in the 
r- and t-directions respectively. The solution was computed both with 
the difference correction neglected and with the correction inserted, in the 
latter case the iterations at each step being stopped when successive 
iterates agreed to six decimal places. As a further check the solution was 
obtained with n — 20, 6f — } and the difference correction incorporated. 
Except in the immediate neighbourhood of the singularity at r Pe 0, 
the solutions agreed to at least three decimal places indicating that the 
values of @(r, t) in Table | are correct to the number of figures given. 

The total time of performing in binary arithmetic the computation 
specified in Table 1 for the case n = 12, 8¢ = 4} was 1 minute without 
application of the difference correction and 2} minutes with the difference 
correction applied. The ACE is a machine having a multiplication time for 
fixed point arithmetic of } a millisecond. 


6. On the adequacy of the finite-difference representation 

In section 3 it was stated that the Crank—Nicolson method is stable. 
This means that, whatever intervals in dt and dr are used in the difference 
equations (8), the rounding errors introduced in the solution at any step 
will not build up in later stages. 





ON A CYLINDRICAL HEAT-CONDUCTION PROBLEM 379 


A related question is the adequacy of the finite-difference representation 
itself. In heat-conduction problems it is not uncommon to have an initial 
discontinuity in function or derivative at a boundary point. Thus in the 
problem of section 2 the value of (20/ér),_, is discontinuous at t = 0. We 
would not expect the finite-difference equations to be meaningful in this 
region for small values of the time and it is often recommended (cf. p. 76 
of (3)) that the difficulty be overcome by using an analytical solution for 
small values of t or by removing the singularity by a transformation of both 
independent variables. 

The numerical results discussed in section 5 indicate that neither of these 
artifices is necessary in the present example if the solution for small values 
of the time is not of particular interest, since the effect of initial discrepan- 
cies in the solution appears to die away for increasing values of the time. 

In this section we examine the validity of this conclusion for the special 
case F(t) = 1 for which an analytic solution of the problem is known 


(ef. (4), p. 176), namely, 
~. 2Jq(7B;) ( en) 
: es | oe oR (19) 
joo (1+ BF). o(B;) ’ Z 
where the 8; are positive roots of 


BJ,(B) = Jo(B) (20) 


arranged in increasing order of magnitude and J, and J, are Bessel functions 
of the first kind. We now compare the solution (19) with the analytic 
solution of the difference equations (8) representing the partial differential 
equation. It is convenient to multiply equation (8) throughout by 2 for 
i = 0. This clearly does not change the solution. 

We consider first the solution of (8) when the difference-correction terms 
C,(@;) and C(¢;) are neglected. In matrix notation equation (8) has then 


the form AO(t | 5t) BO@(t), (21) 
where A and B are square matrices of order (n-+-1) independent of ¢ and 
6(t), 6(¢-+ dt) are column vectors with components 6; and ¢,, respectively. 
Clearly at time m dt we have 

(on 5t) = (A By" 0(0). 
Suppose the matrix A~'B has eigenvalues A, and eigenvectors v;. We may 
express 


6(0) = > CV; (23) 


i 


where the c; are constants. Equation (22) can now be rewritten as 


. 


O(m dt) = 3 CAP V;. 


i=0 
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If A, is the eigenvalue of largest modulus and c, + 0 the right of this 
equation reduces to its first term when m is large. On the other hand, for 
large ¢ the right-hand side of equation (19) also reduces to its first term. 
We therefore have to consider the validity of the equations 

Bi 

Z 


Bit 


Ar exp| Z)" 


l 
that is, —logA, = 
St Ao 


2S, (7, Bo) 
(1 T Bo)o(Bo) 


and Co Vo 
Let us put Z(dr)? 6t — s. Then 
A (C+ 2s]), B = C—2sI, 


where I is the unit matrix of order n+ 1 and C is the square matrix of 
order n+ 1 whose non-zero elements are given by 


r 4 —4 
} 2 -i 


1 ! 
(1 2(n I (I+ 57, 1, 


24-26r(1 + 46r) J 








Note that C is independent of 6¢. The equation defining A, and v, is 

(A B—A,I)v,; = 0 (29) 
or equivalently (B—A, A)v, = 9. (30) 
On using the relations (27) we obtain 


{- 2a(1—A,),) , 
C I'v 0, 31 
we Vi (31) 


1 
showing that the vectors v, are the eigenvectors of C and the A, are related 
to the eigenvalues yx, of C by the equations 
28(1—A,) 28 —p; 


A; 
l + As s+ by 


ey (32) 


The time dependence of the solution enters through equation (32) since 
8 depends on 6. 
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The eigenvalues yp; of C are real and positive and we may suppose that 
By < Byes, fori = 0, 1,..., n. We see then that A; > A;,,. Note that for any 
values of 5¢ and n, we obtain the stability criterion A, < 1. In order that 
A, should be the dominant eigenvalue, as required in setting up equations 

25) and (26), we must ensure that A, > (A, .. This leads to the condition 
. 2Z 
8 > 4, (uou,), Or, equivalently, dt < N ; (33) 
\ (Hon) 

which gives us a restriction on the permissible values of 4¢ for a given value 
of n. In practice the above condition is only restrictive on the initial value 
of dt, for if we choose an initial time interval well inside the critical value 
defined by equation (33) the component of v,, in the solution (24) will be 
reduced after several steps to negligible proportions. We may then increase 
the time interval beyond the critical value. Associated with the new time 
interval will be a different set of eigenvalues, say A,, but the same eigen- 
vectors V,, since in the present problem these are independent of the time 
interval. Thus if we take m steps at an interval 4f followed by p further 

steps at interval 67’ we obtain : 
O(mdt-+pST) = ¥ c,AMAPY,. (34) 

i-0 
If c, A™ is negligible compared with c,Az" it is no longer necessary to have 
\, < Ay. In the subsequent analysis, however, only a fixed time 
interval will be considered. Results for a variable time interval are easily 
deduced from equation (32). 

Tables 24 and 2 relate to the analytic solution (19) with Z = 25 whilst 
Tables 3, 4, and 5 relate to the solution (24) of the difference equations. 
Table 3 gives selected values of «; for n = 10 and n = 20, Table 4 gives 
values of (1/dt)log A, obtained from equation (32) with dt — } and Table 5 
gives the components of ¢) V, and c, Vv, atr = 0,4, lforn = l0andn = 20, 
Table 4 should be compared with Table 2a, and Table 5 with Table 28. 
The condition, specified by equation (33), that A, be the dominant eigen- 
value gives 6¢ < 1-81 for rn = 10 and dt < 0-90 for n = 20, Thus in 
doubling the number of points we have halved the maximum permissible 
value of 5f according to this criterion. Douglas (5) has given a similar result 
for the Crank-Nicolson procedure applied to the equation @6/et = ©*6/ex*. 
He shows that under certain conditions the method converges to the true 
solution as 6t + 0 if the ratio 5t/dz is kept constant. 

When the difference correction is incorporated equation (21) is re- 


placed by (A+D)@(t-+5t) — (B—D)@(t), (35) 


where D is a matrix arising from the difference correction. For definiteness 
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TABLE 2a TABLE 2B 








2] o(rBe) 2/,(rBo) 
(1+ 83) Jo(By) | (1+ BF) /o(B,) 
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TABLE 3 TABLE 4 








} 1 | 1 
| log Aj (nm 10) Fe bog Ay (n 20) 
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O'0H308 1 0706308 1 


066146 s o'66577 s 























we shall assume that fifth and higher-order differences are neglected in the 
formation of ((4,). The computation of C,(6,;) for ¢ 0, 1, n—l,andn 
requires a knowledge of @_,, 6_,, 4, .,, and @, ., which he outside the normal 


region of integration. By symmetry @_, = 6, 0_, = 6,; also we suppose 
that @,., and @, ., are computed from the extrapolation formula 


6,45 60,— 150, ,+- 208; »— 156; 44 66,4 0, 5 (a n,n+1) (36) 


obtained from the condition that V°,_, and V0, should vanish. Then the 
non-zero elements of 12D are given by 


16 4 


2 
2 6 
j 


3 6 
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The analysis leading to equation (32) can be extended with C now replaced 
by C—D and Tables 6-8 correspond to Tables 3-5. To ensure that A, is the 


dominant eigenvalue we now require dt < 1-59 for n = 10. 


TABLE 6 TABLE 7 TABLE 8 





4 
0° 28965 
005834 
+ O11 327 


py in 10) “oVe 





- log Aj (nt = 10) 
C’O1S77 0 1*20709 0 

rrogior 6 
o-7 7609 6 


o16625 8 o'063708 1 


H2411S 3 066656 9 




















We see from these results that equations (25) and (26) are asymptotically 
satisfied for large t,so that eventually the solution of the difference equation 
must closely match that of the differential equation. The initial dis- 
| has no lasting effect. We further notice that the results 
10 with the difference correction incorporated are better 
approximations than those for n = 20 without the difference correction. 
On the other hand the time taken to solve the equations for n = 10 with the 
difference correction inserted is as great as the time taken to solve those for 
n — 20 with the difference correction omitted. This comparison, however, 
depends on the linear nature of the problem and the situation could well be 
different in 


continuity at r 
obtained for n 


a non-linear case. 


TABLE 9 


n 10, 


bt = 3 





In flerence correction omitted 


Difference correction inserted 





o*e(0, f) 


10%e( 4, 0) 


10%e(1, f) 


10%¢e(0, f) 


10%}, 


10%e(1, f) 





301 
665 
490 
275 
115 

9§ 
$7 
14 

4 





rit 
55 
13 
21 
14 
13 
6 
° 
° 


743 
124 
s6 
151 
204 
225 
126 
37 
11 


108 
93 
136 





1,673 
109 
66 

4° 

26 


8 
6 
3 
1 

















Table 9 illustrates the behaviour of the difference «(r,t) between the 
analytic solution and the numerical solution for the case n = 10, 5t = }. 
We see that ¢(r, t) becomes substantially smaller for large values of the time. 
Of course, since the steady-state solution is independent of the initial 
temperature distribution any process which tends to the steady-state 
solution will appear to become increasingly accurate to fixed-decimal 
working if we consider sufficiently large values of the time. However, 
Table 9 shows that the numerical solution with the difference correction 
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inserted is in error by less than six units in the fifth decimal place for t > 5. 
Atr = 0,t = 5the analytic solution is 0-870174 which is still far removed 
from the steady-state value of zero. 


We have verified that the solutions obtained numerically by the pro- 


cedure outlined in section 4 are in agreement with the analytical solution of 
the difference equations given by (24). 


7. Conclusions 

Problems of the type specified in section 2 are eminently suitable for 
direct numerical solution on automatic computers by use of the Crank-— 
Nicolson implicit method. Although the finite-difference representation 
is inadequate in the neighbourhood of the singularity at the boundary the 
errors arising from this source in the present problem diminish with in- 
creasing time and no special starting procedure is necessary to deal with the 
singularity. The accuracy of a given solution can be checked by recomputa 
tion using different intervals in the variables. 
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